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The generalized density matrix (GDM) method is used to calculate microscopically the parameters 
of the collective Hamiltonian. Higher order anharmonicities are obtained consistently with the lowest 
order results, the mean field [Hartree-Fock-Bogoliubov (HFB) equation] and the harmonic potential 
[quasiparticle random phase approximation (QRPA)]. The method is applied to soft spherical nuclei, 
where the anharmonicities are essential for restoring the stability of the system, as the harmonic 
potential becomes small or negative. The approach is tested in three models of increasing complexity: 
the Lipkin model, model with factorizable forces, and the quadrupole plus pairing model. 
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I. INTRODUCTION 

A long-standing question of microscopic description of 
nuclear collective motion belongs to the class of problems 
which are left behind by the advancing army that cur- 
rently is mostly interested in new frontiers, in our case, 
in drip line physics. Meanwhile, we still lack a system- 
atic theory based on first principles and inter-nucleon 
interactions that would allow us to fully understand nu- 
merous collective phenomena in the low-energy region of 
medium and heavy nuclei and satisfactorily describe the 
data. In relatively light nuclei, the shell model (what is 
nowadays called configuration interaction) with effective 
nucleon-nucleon forces usually works well although even 
here the abundant numerical results sometimes require 
some kind of model interpretation. In heavier nuclei, the 
necessary orbital space is too large for direct numerical 
diagonalization. 

Phenomenological models frequently work well, first of 
all the geometric Bohr Hamiltonian [l|, y] and the inter- 
acting boson model (IBM) Q- However, the relation be- 
tween their parameters and the underlying microscopic 
structure remains uncertain. Moreover, some assump- 
tions of such models turn out to be unreliable. For ex- 
ample, the identification in the IBM of the prescribed bo- 
son number with the number of valence fermionic pairs 
breaks down in the attempt to explain very long "quasivi- 
brational" bands extended, without considerable changes 
in spacing, up to spin values much greater than the finite 
boson number would allow, see for example the ground 
state band in ^^'^Cd close to the equidistant ladder up to 
,r = 28+. 

The microscopic theory is relatively successful in well 
deformed nuclei. Various mean-field methods, includ- 
ing the modern energy density functional approach [J, Q 
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with pairing, indicate regions of nuclei with clearly pro- 
nounced deformed energy minima. With the microscopic 
definition of shape, one can calculate the moment of iner- 
tia by the cranking model and the generator coordinate 
method, construct rotational bands built on different in- 
trinsic configurations and explain back-bending and sim- 
ilar phenomena [6|. 

In our opinion, the status of microscopic theory is still 
underdeveloped with respect to spherical nuclei, espe- 
cially in the case of the presence of a low-lying collec- 
tive mode. The standard way of defining such modes is 
based on the quasiparticle random phase approximation 
(QRPA). This is essentially the harmonic approxima- 
tion that determines the frequency and two-quasiparticle 
structure of the collective phonons. If the multipole cou- 
pling is strong, the collective mode has a large ampli- 
tude, the frequency falls down, and the QRPA reveals 
the instability. In reality, this is not necessarily a point 
of phase transition. Rather, this is the region of strong 
anharmonicities outside of the reach of the QRPA. Phe- 
nomenologically, this can be described by a special choice 
of potential and rotational parameters in the Bohr Hamil- 
tonian which are close to the 0(6) limit of the IBM with 
a gamma- unstable potential. Currently we do not have 
a reliable microscopic approach to quantify collective be- 
havior of this type. Another practically important ques- 
tion related to anharmonicities is the mode-mode cou- 
pling. The coexistence and interaction of soft quadrupole 
and octupole modes are relevant, for example, to the 
search of mechanisms for many-body enhancement of the 
nuclear Schiff moment and the atomic electric dipole mo- 
ment 0]. 

Instead of the direct diagonalization of the primary 
nucleon Hamiltonian, it seems reasonable to work out 
a procedure for the microscopic derivation of the effec- 
tive collective Hamiltonian. Typical collective states can 
usually be identified by their quantum numbers, low en- 
ergies and large transition probabilities. Being intercon- 
nected by large matrix elements of corresponding col- 
lective operators they form a collective subspace of the 



total Hilbert space of the system. In the case of a soft 
muhipole mode, it is often possible to label the empiri- 
cal levels by the phonon quantum numbers, even if their 
energies and transition rates noticeably differ from the 
predictions of the harmonic approximation. This differ- 
ence results from anharmonic effects which still keep the 
geometric nature of the mode. Therefore our approach 
will be to develop the road to a consistent mapping of 
the underlying nucleonic dynamics onto that inside the 
collective subspace. 

The idea of this approach goes back to the boson ex- 
pansion technique suggested long ago Q; a detailed re- 
view of work in this direction can be found in [9| . The 
formalism of the generalized density matrix (GDM) re- 
formulating earlier work [l^ by Kerman and Klein seems 
to be the most appropriate for our goal [11 -[l4l. T his for- 
malism was applied to collective rotation lll.Il3l . [isl [la | 
and large amplitude collective motion [17h20|| generaliz- 
ing the time-dependent mean- field method [2l| . Here we 
apply the GDM approach to collective vibrations in soft 
spherical nuclei. 

The generalized density matrix R12 = a^ai is the set 
of operators defined originally in the entire Hilbert space 
[1 and 2 here represent a complete set of single-particle 
(s.p.) quantum numbers]. The microscopic Hamilto- 
nian provides exact operator equations of motion (e.o.m.) 
for this set. Taking matrix elements of these equations 
between the states of the collective family we map the 
equations onto the dynamics of the collective operators 
inside this family. The choice of the collective Hamilto- 
nian should be quite general dictated by the type and 
symmetries of collective motion under study. Compari- 
son with microscopic dynamics determines the collective 
parameters. The lowest orders give naturally the mean 
field [Hartree-Fock-Bogoliubov (HFB) equation] and the 
harmonic part (QRPA). Next orders determine anhar- 
monicities. These higher order terms are not assumed to 
be perturbative, they are separated only by their oper- 
ator structure in the collective space. Simple estimates 
|23 . [23j show that in many generic cases the quartic an- 
harmonicity with respect to the quadrupole coordinate 
plays an important role. In fact, this was earlier con- 
firmed by specific realistic applications '2,4^ of the phe- 
nomenological anharmonic Hamiltonian; ^°"Pd is proba- 
bly the clean example of such dynamics. 

We start with the discussion. Sec. [Hi of the general 
procedure of the GDM method. In Sec. IIIII we consider 
systems near the critical point (small RPA frequency oj^). 
Sec. HV] and Sec. |V] are devoted to the Lipkin model 
and factorizable force model, respectively, which tradi- 
tionally serve as a testing ground for various theoretical 
approaches. Sec. IVII discusses the GDM method applied 
to realistic nuclei with pairing and rotational symmetry. 
In Sec. IVIII we give the results for a quadrupole plus 
pairing Hamiltonian, with a semi-realistic numerical ex- 
ample. Sec. IVIIII summarizes our method and discusses 
future working directions. The details of calculations are 
given in the Appendices. 



II. THE GENERALIZED DENSITY MATRIX 
METHOD 

In this section we reveal the essence of the GDM 
method, in a simple system without complications due to 
rotational symmetry and pairing correlations. A single 
collective mode is assumed; the case of multiple modes 
is discussed briefly in Appendix |L] The main result, be- 
yond the well known HF equation and RPA, is a relation 
(1441) involving cubic and quartic anharmonicities. 



A. Preparation 

The starting point is the effective microscopic 
fermionic Hamiltonian 



-ff = ^ Zi2a\a2 + T X! Vi234a\ala3a4. 
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We find it convenient for H, Z12 and V1234 to be dimen- 
sionless; in other words H is measured in some unit of 
energy. We have assumed in eq. ([Ij a two-body force, in- 
clusion of three-body forces is discussed in Appendix [A] 
In accordance with the discussion in Sec. U we assume 
that H has a band of collective states {jCi)} character- 
ized by low energies and large transition amplitudes. We 
assume that there exists a reference state |$), a collec- 
tive mode operator A^ = {a — iTr)/^/2 (a, n are collective 
coordinate and momentum), such that approximately 

[a, tt] = i, (2) 
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Eq. ([2]) says that A^ is effectively a boson operator. Eq. 
([3]) says that the collective band {jCi)} can be built by 
repeated action of A'' or A on the reference state ]$). 
Later ]$) will be identified as the HF ground state. Eq. 
(j4]) says that within the band, the effect of the fermionic 
Hamiltonian can be approximated by an expansion over 
the bosonic operators, where we keep all time-even terms 
up to quartic anharmonicities {a is time-even, tt is time- 
odd). 

Now our goal is to map the exact e.o.m. in the full 
Hilbert space onto collective dynamics inside the band 
subspace. We will use contractions and normal ordering 
of operators. They are defined as: 



A'B' = ($1^5]$) 
N{AB\ =AB- {<i>\AB\<P) 



(5) 
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Without paring, the reference state |$) has a definite 
particle number, 

ma\a2\<i>) = P21, ($|aia2|$) = ($|al4l*>=0. (7) 

p is the usual single-particle density matrix. Normal or- 
dering of more than two operators is defined by the Wick 
theorem: 

a| 030304 = iV[oi 030304] -I- p4iN[ala3] — p3iiV[o204] 

-p42iV[o|o3] + P32N[a\a4] + P41P32 - P3lP42- (8) 

Equivalently, normal ordering puts quasiparticle creation 
operators to the left of annihilation operators. 

The generalized density matrix operator is defined in 
the full space as 



i?i2 = otfli ^ P12 + N[alai] = P12 + R.^2- (9) 

The Hamiltonian ([T]) in the normal ordering form is 

H = mH\^) +J2fi2N[a\a2] 



t.t. 
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where we have introduced the self-consistent field opera- 
tor 

W{R}i2 = J2 ^1432i?34, /l2 - ^12 + W{p}i2, (H) 
34 

and mH\<^) = ^12(^12 + \W{p)i2)p2i is the average 
energy of the reference state. 

The exact e.o.m. for the density matrix operator in 
the full many-body Hilbert space is 

[Ri2,H]^[a\ax,H]^ 

[f,p]l2 + [f.R'']l2 + [W{R''},p],2 

+ 2 Zl(^i345^[4a3a4a5] - V5432iV[o^ola3ai]). (12) 

345 

Since we are only interested in the band subspace, we 
take matrix elements of eq. (IT^ between two collective 
states: 

(a|[i?i2,iI]|C,) = 
m[f,p\l2 + [I,R'']l2 + [W{R''},p\,2 

+ l{^{yiM^N[alala4a^]-V^432N[ala\a3a{\) |Q)(13) 
345 

We assume that within the band the effect of i?i2 can be 
approximated by a boson expansion: 



{Ci\R^2\C.j) = (Ql P12 -f ri^2 ^a + rl^2 ^TT 



(10), 
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where we keep explicitly terms up to quartic anharmonic- 
ities. A convenient normalization is: a term with m of a 
and n of TT has a factor of 1/ [mn); each anti-commutator 
gives an additional 1/2. Similarly, for A^[a4oJo20i] we 
have 



xloto20i]|Q) = (CI \r'^la^ + \r^l^^ 

+ 2 12341"^' ^J + 3^1234*^ + 2^^1234- 



-\rg'lW^n} + \r[ll{a,n^} + ... |Q) 



where we have assumed that the expansion starts from 
a^, TT^, {a,7r}, as explained in Appendix [Cl Now the 



r.h.s. of eq. (jl3|) is written as an expansion over boson 
operators. 

The l.h.s. of eq. (|13p is approximately given by: 
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12 



.(10), 



{Ci\[Ri2,H]\Cj 



Eq + — a' + ..]\C. 



3h 



where we have restricted the intermediate states (be- 
tween i?i2 and H) by those of the collective subspace 
{|Ci)}, since H is a, collective operator: the matrix ele- 
ments of H connecting the collective band with states of 
a different nature are small. This is the main approxima- 
tion of the method; influence of the neglected "environ- 
ment" states can be later accounted for with the use of 
statistical assumptions [2^. After calculating commuta- 
tors like [tt, a^] = —2ia, the l.h.s. is written as a boson 
operator expansion. Then we equate in eq. (J13l) l.h.s. 
and r.h.s. coefficients of the same phonon structure: 1, 
a, TT, a2^2 ... The resultant equations are examined be- 
low. 



B. Zero Order: Mean Field (Hartree-Fock) 



Terms without a or tt in eq. (J13p give 

[/,P]12=0. 



(16) 



Thus / and p can be diagonalized simultaneously in some 
s.p. basis: 



/12 = (^1261, P12 = Si2ni, 



(17) 



providing mean-field s.p. energies and occupation num- 
bers. We will always use this s.p. basis. If we restrict 
the reference state |$) to be a Slater determinant, then 
the occupation numbers ni can be only or 1; in this 
case eq. ([16]) is the usual HF equation, |$) is the HF 
ground state. More general choices, such as the thermal 
ensemble, are also possible. For future convenience we 
define 



' 12 



612 = 61-62, ni2 = ni-n2. 



(18) 



We assume that degenerate s.p. levels have the same 
occupancies, 



ei = 62 ^ ni = n2, 
but the reverse is not necessarily true. 



(19) 



C. First Order: Random Phase Approximation 

Terms linear in a and tt in eq. (J13p give 

tt: zr(io) = [/,r(oi)] + [u;(oi),p], (20) 

a : - za;2r(°i) = [/, r^io)] + [w;(i") , p] , (21) 

where w^^o) = M^jr^io)}, and u;(°i) = M^jr^oi)} ^^^.g ^^le 
corresponding components of the mean field. This is the 
set of RPA equations. The formal solution is 



„(io) 



ni2 



-zuj^ w^^^'^ + ei2w[l% (22) 

(01) -_ ^^12 L- „„(10) . ^ ,„(01)i 



'' "(ei2)2-c.2 
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f, ^, [.^r+ei2^r]- (23) 

^ei2j — w 



Note that r^i") and r^o^' have only ni ^ n2 matrix ele- 
ments. From eqs. (dJ), (HH and ([25)1 we obtain a linear 
homogenous set of equations for w^^") and w^"!); 



(10) 
"34 



-E^3 



?^12 



3214 



12 



(ei2)2 - W 



W. 



(01) 
34 



E^= 



»^12 



3214 



(ei2)2 -w- 



[-*a.24°/)+ei2u;if],(24) 
:Mf +ei2«;l°2'^].(25) 



Introduce the matrix M: 



M = 






(26) 



where (ea ^64 , ei 7^ 62) 



^("34), (12) = ^("34), (12) = <5(12),(34) - ^32147^3^^^,(27) 



(ei2)2-c.2 



^(34), (12) — ^ ■ ^(34), (12) 



^3214 .^.2"" 2 ^(28) 



(612)^ -CJ 



in which (5(i2).(34) = '5i3<524- Then the 63 ^ 64 part of eqs. 
dMl) and d^ is written as 



M ■ 



w 



(10) 



0. 



(29) 



^ . u;(01) 

Non-zero solution requires a zero determinant: 

Det[M] = 0. (30) 



Eq. (pO|) is the RPA secular equation determining the 
harmonic frequency w^. 

By eq. pop the transpose matrix of M, Af^, has a 
zero eigenvalue. Assume the corresponding eigenvector 

is (A34,X34): 



M 



T I ^\ _Q /A 



In other words, the row vectors of AI are linearly depen- 
dent. A and x are used later. 

The normalization of A^^\ r(°i) can be fixed by the 
so-called saturation principle as explained in Appendix 
11 



? [(ei2)^ 



-n-12 



,212 



/• (10) 
■(iwl2 



(01)w 2 (01) 
612^2 )(^ ^"21 



I ■ (10)^ 

+ 4621^21 )• 



(32) 



Starting from the next (second) order, ri^^l of eq. (fT5)) 
begin to appear in the e.o.m., the saturation principle is 
also used to express them in terms of r^™"' of eq. (|14p . 
as explained in Appendix IB] 



D. Second Order: Cubic Anharmonicity 

Using eqs. (|B1HIB13| . the second order terms, a^/2, 
7rV2, {a,7r}/2, in eq. ^ give 

- 2iwV(") - 2iA(30)r("i) = 
[/, r(20)] + [w;(20) , p] + 2[ii;(lO) , r^")] , (33) 
-i^V("2) + zr(20) + zA(i2)r(io) = 
[/,r(")] -H M"),p] + Mio),r("i)] -I- [ii;(oi),r(i")], (34) 

2ir(") - iA(i2)r(oi) = 
[/,r(02)] + [u;(02)^^] + 2[u;(oi),r(oi)]. (35) 

First we consider the ei = ei' matrix elements rj^-^, , 



r!^;^, and rl^' . As one can check, eqs. ([55)1 and ([55)) 
give the same r]^;^, . But r\^, and rj^;^, are not fully 



determined since eq. ([34|) determines only the difference 

'^^^ii' ~ '"11' ■ ^^ ^-'^ them by the saturation principle, 
the ei = 61/ matrix elements of eqs. (jB6p and (jBSP : 



rl7=-z[p,r(io)],,„ <f^=z[a;,r(oi)],,„ 

= -z[p,r(oi)]ii,. (36) 



r(")=zfx,r(io)l 



It is straightforward to show that eq. ([5S)) is consistent 
with eqs. §MB- 

Next the ei ^ 62 matrix elements r^2 : ^12 and rj^j 
are solved formally in terms of u;(20/02/ii) ^^^^^ ^(30/12) 
from eqs. p3ll35p . Then by eq. pT|) we obtain the Hnear 
set of equations for w^^^\ w(°^) and w^^^\ 



E. Third Order: Quartic Anharmonicity 



^ (A^ X^)-M = 0.(31) 



Using eqs. (jB14IIB17|) . the third order terms, a'^/3, 
{a2,7r}/4, {a,7r2}/4, n^/3, in eq. dH]) give 

_ ^^2^(21) _ 3,A(30)r(") - 3iA(«)r(oi) ^ 
2 

[/,r(30)] + [u;(30),p] + ^[^(20)^^(10)] + ^^^(10)^ ^(20)]^ (37) 



+iA(22)^(10) ^ [^^^(21)] ^ [^(21)^^] + [^(20)^^(01)] 

+ [«;(oi),r(20)]+2[ii;("),r(i")]+2[«;(io),r(")],(38) 



[/,r(12)] + [^(12),p] + [^(10)^^(02)] 

+ [^(02)^^(10)] ^2M"),r(°i)]+2[«;(oi),r(")], (39) 



^^(12) _ ^ A (12) (02) , 3jA(04) (10) 
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[;,,(03)] 



[^(03)^p] + 3[^(02)^^(01)]^3^^ 



(01) ^(02)1 



(40) 



The ei = ey matrix elements rly can now 

be found in terms of the fower order quantities from 

(30/21/12/03) 



eqs. p7M0p . The ei ^ 62 matrix elements rjj 
can be calculated formally in terms of u''^^^/^^/^^/'^'^^ and 
^(40/22/04) fj.Qj^ gqg ([3711401) . By eq. ^ we obtain a 
linear set of equations for w^'^"', w'^^^ w^^^^ and z/;^°^'. 



However, this set is not linearly independent. Thus 
we have a solvability condition. To see this, keeping 
only (r/u;)(30/2i/i2/03) and A(40/22/04) tg^jns^ eq^ ^ 

+ iw^xeq. dSni) gives 

_ jC^2j 1^(21) ^ ^2^(03)] ^ [^^ ^(30) _^ 1^2^(12)] 



,(30) + 1^2^(12)^^] ^3^^(40)^(01) 
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^xeq. 



w xeq. (|40|) gives 
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The variable parts of eqs. (HIT) and (jH]) have the same 
structure as the RPA equations (|20p and (HH). Introduc- 
ing temporarily x — r^'^o^ + iw^r^^^) , y = ir'^^^+w^r^^^^ 
we can solve a;, y in terms of VK{a:}, VK{y}. Using eq. 
(fTTj) to obtain linear equations for T4^{x}, W^{y}, the 
63 7^ 64 part is written as 



M 



W{x} \ _ f A 
* W{y} ^[b 



(43) 



where the matrix M is defined in eqs. (|26II28P : A, B 
consist of A^""*), A^^^^ and A^^"*) terms. Multiplying eq. 
(|43)) from left by (A-^ x^) and using eq. ([3T|) we come to 
the solvability condition: 



3A(«)^ E[X34-ei2A3 
12 e37^e4 



JOl) 



14: 



(ei2)2 



+ ^A(-)EEF3214 
12 637^64 

2\ r (10) • (oi)i r (10) ■ 2 {oi)-| 



(ei2)2 



„(io) 
'12 



+3C.'A(0") V V [a;2A34 - ei2X34] • V3214^^^ J = . . .(44) 

l^et^e, (ei2)2-a;2 

The r.h.s. ". . ." contains only the lower order quantities, 
including A^'^"^ and A'^^^^ On the l.h.s. the coefficients 
of the A(40), A(22)^ A(04) terms are of order 1, u^, uj'^, 
respectively, although this might not be obvious from eq. 
(HH). It follows from examining the expressions of Af , A 
and B in eq. (|43p. This point will be important for the 
discussion in Sec. IIIIl 

At the current stage, the cubic and quartic anhar- 
monicities are not completely fixed, we find only one re- 
lation (im constraining them. However, we are able to 
obtain A'-'^"-' and A'""^' near the critical point a;^ « 0, with 
important applications, as will be explained in Sec. IIIIl 
Even with this limitation, eq. (|44|) is useful. One could 
fit the ratios of A^™"^ with the experimental data, then 
use eq. (|44p to determine their magnitudes. This is espe- 
cially interesting for the cases with certain symmetries, 
where the ratios are known. Results in this direction will 
be discussed elsewhere. 



F. Self-consistent Hamiltonian Conditions 

If the approach is self-consistent, substituting the so- 
lutions of eqs. (jM)) and (fTS]) into eq. ([T]) should provide 
eq. (|3]). Namely, 

Eo = mm^) = E(Zi2 + ^w^{p}i2)P2i, (45) 

12 
A (mn) \ ^ Imn) t ^— \ (mn) /ap\ 

A^ ^=2^eirii ' + 72^1^234^32/. (46) 

1 1234 

We checked eq. (|46l) explicitly up to the cubic an- 
harmonicities. Concerning the quartic anharmonicity, 
we checked the combinations A^"'"-' 



1, ,2 
2 



cjj-AC^^) and 



[^,2^(04) _ 1^(22). since from the fourth order e.o.m. (not 
listed) only r[^ 



(40) 



iw^rii ' and u'^r) ' 



'11 



'11 



1„(22) 
2'll 



are 



determined, similar to the situation in the cubic order. 

In summary, this section discusses the general proce- 
dure of the GDM method. The exact e.o.m. for the 
density matrix operators are mapped onto the collective 
subspace by taking matrix elements between states of 
this family. Comparing terms with the same phonon op- 
erator structure, order by order, we get equations for the 
GDM. In each order, the GDM is solved from a set of 
coupled linear equations in terms of lower order quanti- 
ties. The bosonic Hamiltonian coefficients A^™") appear 
as parameters in the solution. 

At the current stage the anharmonicities are not com- 
pletely fixed; we find only one relation (|44p involving cu- 



bic and quartic anharmonicities, appearing in the third 
order as a solvabihty condition. In the next section, we 
will show that the cubic potential A^'^"-' and quartic po- 
tential A'^'*°^ can be determined in a special case - around 
the critical point w^ « 0. 



III. SYSTEMS NEAR THE CRITICAL POINT 

Anharmonicities become important when the har- 
monic potential ixP'c? jl becomes small or negative. This 
is the case in many realistic medium and heavy nuclei 
away from magic numbers [26| . The quartic potential 
AC**") and higher terms restore the stability of the system. 
At the same time, the system can be deformed by odd 
anharmonicities; the potential is flat at the bottom, or 
gamma- unstable. Near the critical point uP' ~ 0, we are 
able to determine the cubic potential term A'^^") and the 
quartic potential term A^'"'^. Deformation due to A^"^"^ 
will be studied separately. In this work we concentrate 
on the case of small A^^°) , consistent with the idea of soft 
spherical nuclei. 

We make an assumption in the spirit of Landau phase 
transition theory: in eq. Q, the leading potential term 
w^a^/2 vanishes at the critical point, while other higher 
order terms A^^™") remain finite. Taylor expanding A*^™") 
over w^, 



^(mn) ^ p^(rnn) ^ ^i^^n) ^2 ^ ^^nn) ^4 ^ 



(47) 



, (rnn) 



the leading constant term Ac is finite. 

Near the critical point the stability of the system is 
restored by higher order anharmonicities, e.g. A^'^^^a^/i. 
Thus {Ci\a\Cj), {Ci\n\Cj) ... are finite. Consequently 

r^™" in eq. (|T4|) is finite, since the l.h.s. {Ci\alai\Cj) is 
finite. Again we call the finite leading constant term in 

a iaylor expansion r^-^^2 ■ 

We can obtain Ac from eq. (|44)) by keeping only 
leading constant terms (neglecting terms with w^, w'* ...), 
as explained below eq. (|44p . Another approach is pos- 
sible: neglecting w^ terms earlier, in each e.o.m. For 
convenience, we use = instead of = if an equation is cor- 
rect in constant terms but not in oj^ terms or higher. In 

this way we determine Ac , eq. (|53p. as well as Ac , 
eq. dSSl). 



A. RPA 

Keeping only the constant terms of eq. (I24[) we have 

1^3214 Wi2 ' . (48) 



W 



(10) ^\^T/ "12 ^JIO) 



34 



12 



ei2 



Defining a square matrix 



n-12 



-0(34), (12) = ^(12), (34) - "1^3214 , (63 ^ 64, 61 7^ e2),(49) 

C12 



the 63 7^ 64 part of eq. (|48p is written as Dw^^'^^ = 

0. Since the quantities w\i2 ^^ not vanish, we have 
Det[D] = 0. Thus D^ ^ the transpose matrix of 13, has 
a eigenvalue. More accurately, ti^ has an eigenvalue 
of order w^; because Det[D], the product of all eigenval- 
ues of D"^, is of order us^ . Assume that the eigenvector 
corresponding to this eigenvalue is 7734: 



D^r\ = 77 = ^ jfD 



(50) 



B. Cubic Anharmonicity 

Keeping only the constant terms of eq. ([33]) . 
- 2zA(3°)r(°^) = ei2r(f - ni2^1f + 2[u;(io), r(i°)]i2, (51) 
and calculating W34 from eq. (|5ip . the 63 ^ 64 part is 



written as 



D ■ «;(20) ^ c, 



(52) 



where D is defined in eq. (P^)) . C contains A'^'^") and 
lower order quantities. Multiplying eq. (1521) from left by 
rf^ and using eq. (|50p we obtain 



A(30) . ^ J2 VUVS2U^ 
ei#e2 63^64 



.(01) 



ei2 



i ^ ^ ?734^3214- 



zi;(io),r(io)]i2 



ei#e2 635^64 
1 
'2 



ei2 



I E E m.V,ru[pM'%,. 



(53) 



ei=e-i> 637^64 



where p is given in eq. (jB7p . Eq. (|53|) gives Ac +0{uj^). 
Then w\2 is solved from eq. ([5^ with an overall factor 
still undetermined. 

Similarly, from eq. p4p we obtain an equation D ■ 
w^^^> = . . . Multiplying it by rj^ we obtain 



A^^'^ E E ^34^3214^ 



.(10) 



£15^62 635^64 



ei2 



E>r^ ^. "-12 (20) 

Z^ %4V^3214 7^7^U'l2 



£15^62 635^64 



'(ei2)^ 



2zA(3°) J2 E ^34^^3214 



JOl) 

'12 

(ei2)2 



w;(i"),r(io)]i2 



r, 

'« E E '^34^3214- 



£17^62 637^64 

+2 E E %4 1^3214- , „ 
^ ^ (612)-^ 

61/62 637^64 

y;(io)^^(oi)]^2 + M°i\r-(io)]i2 



617^62 63/64 



ei2 



- E E %4y3i'i4[^(''^\r(°i)]ii'.(54) 
61=61/ 63/64 



(20) 



From eq. (j54[) the undetermined overall factor in wl 
is fixed as a function of A'^^^ Then from the equation 
D ■ w^^^^ = ... we solve for w-fj as a function of A'-'^^\ 
with an overall factor still undetermined. After doing 
similar manipulation on eq. psp . the undetermined over- 



(02) 

Wn, IS 



all factor in wjj is fixed as a function of A*^^^) 

solved as a function of A^^^^^ with an overall factor still 

undetermined. 

In summary, there remain two undetermined param- 
eters in this order: A'^^^ and an overall factor in uij^g ■ 
We will see them explicitly in the factorizable force model 
(Sec. El). 



C. Quartic Anharmonicity 

Similarly, we obtain from eqs. ([57)1 and ((55)) : 



A(«') . Y: E ''341^321 



JOl) 



ei7^e2 esT^ei 



ei2 



Aa(-) E E m.Vs.^.r^ 



(20) 



= 1=^1/ 637^64 



-A(^») E E ''«4^3214- 



(11) 
12 






^iA(-) E E m.v,,,,r[f 

ei=ej/ <:37^<:4 
r„(20)^^(10),^^r (10)^^(20),^ 



ei2 



(20) (01)i I r (01) (20)i 



+ 2 E E '?34V'3214 
617^62 637^64 

"6 ^ ^ '734 1^31' 14 ([' 

-^ E E ^34T/3i'i4(K<"\r("')]ii' + [»<^"\r'^^>]ii0.(55) 

ei=ei/ 637^64 

Eq. ([55)1 gives A^ + O(w^). There is one unknown 
parameter A^-'^^^; quantities (r/ui)^^"^ and (r/u;)^^^' 
depend implicitly on A*^^^^. 

In summary, this section fixes the cubic potential A''^'^-* 
([53)1 and the quartic potential A'-'*"-' ([55)1 near the criti- 
cal point uj"^ w 0, by considering the leading terms of 
the e.o.m. Deformation due to A'-'^"-' will be studied else- 
where. Near the critical point, the stability of the system 
should be restored by the quartic potential A^'"'^ if it is 
positive and large. In the following we test this idea in 
three models of increasing complexity: the Lipkin model 
(Sec. HV)) . model with factorizable forces (Sec. El, and 



the quadrupole plus pairing model (Sec. I VII 



IV. LIPKIN MODEL 

We test the GDM method in the Lipkin model [13] 
where the analytical solution is available. As we will see, 
the agreement is perfect (Sec. IIVCI) . Then we discuss 



some problems inherent to the bosonic approach itself 
(Sec. HVDJ) . 



Exact Solution 



In this model, there are two s.p. levels with energies 
±\ (the spacing is the energy unit), each with degeneracy 
57-1-1. The model Hamiltonian contains only "vertical" 
transitions (cr = ±1; Z = 1, 2, ..., J7 -|- 1): 



^ = E 2"^.i"'^'' ^ 2 ^ aLaI,i'a-'^.''«-'^.i- (56) 
The quasi-spin operators, 

J+ = j1_ = Jx+ iJy = E "+i,i"- 



-i,ii 



1 V- 



2^i^f^V/a'^,i' (57) 



satisfy the angular momentum algebra. Using eq. ((57l 
the Hamiltonian (1561) is written as 



(58) 



H = J, + ]^k{JI + J^_) 



and the total quasi-spin J is a good quantum number. 
With the Holstein-Primakoff transformation (HPT), 



J+ = j1 = ^V2J-At^, J^ = -J + A^A, (59) 

where A'^ and A are bosonic creation and annihilation 
operators with commutation relation [j4,A^] = 1, the 
Hamiltonian ([55)) is written as an expansion over A^ and 
A; or a and tt by the canonical transformation 

A — —j={iua + vtt), A^ = —^{—iua + vir), uv = —1. (60) 
V 2 V 2 

Assuming J 3> 1, we keep only the leading order in 
1/ J. Under the choice 



u ■ 



Vl + 2kJ, V 



1 



(61) 



the Hamiltonian becomes 

2 



«.^„.,l.»,^„.,A!p,.. „2, 



with 



1-Ak^J\ A(40)«ku4, A("4)«_k„^ (63) 



Other A*^™") vanishes in their leading order of 1/ J. 
Around the critical point cj^ « 0, 

K« , U « V 2, w~ =, 

2J' ' 72' 

A(40)_1 A^O'^' ~ ^ 

J' 8J- 



(64) 
(65) 
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The GDM Method 



Applying the GDM method to the Hamiltonian ([& 
we have solved for rj^™" explicitly in terms of A*^™"^ fol- 
lowing Sec. ini Below we summarize the main results. In 
the mean-field order, the HF s.p. levels are the same as 
the original s.p. levels. Introducing n = 7ia=-i — ?Ict=i > 
0, where ria are occupation numbers of s.p. levels, in the 
harmonic order the RPA secular equation ([5(11) becomes 



(66) 



cj2 = 1 - {nnVtf 



In the quartic order, the solvability condition (|44|) be- 
comes 

3 A(4") + A(22) + 3c^2A(04) ^ 
^ _L (A(30))2 + 2^ A(30)A(i2) + ^(A(i2))2 



12 



l-w" 



n{a + l) 



(67) 



C. Comparison with Exact Solution 



The quantum number J is found from eq. ([571 



J=|J. 



r^ + i 



|ni 



n_i| = 



n{n + \) 



(68) 



We assume 2J = n(f7 + l)^l. In the harmonic 
order, the RPA secular equation (1551) agrees with the 
HPT frequency equation (|63p. In the quartic order, the 
HPT solutions ([63]) satisfy our solvability condition (l67l) . 

The diverging behavior of eq. (|67p around the critical 

point w^ ~ gives Ac . The Ac term must vanish as 
seen from the presence of the term -y -^ (A^'^"^)^, which 
is the only one divergent as o;"^. Equating the l.h.s. and 
r.h.s. diverging terms ^ uj~'^ we obtain 



Al^o) 



n{Q. + l) 



(69) 



This agrees with the HPT solution ((65l), A^^^) 



If we follow the procedure in Sec. IIIIl we obtain 



n(n+i 

the same result ([M)) 



D. Numerical Diagonalization and Discussion 

Here we discuss some problems inherent to the bosonic 
approach itself. The bosonic Hamiltonian ((4]) is usually 
diagonalized in the infinite phonon space; practically the 
space is enlarged until convergence is reached. However, 
there exists a maximal phonon number, close to the ac- 
tive valence particle number in the system. Applying 
the phonon creation operator A^ too many times to the 
ground state, we run out of valence particles. We will 



call this finite phonon space "physical space". Only if 
e.g. the first excitation energy has reached convergence 
within the physical space, it is valid to formally enlarge 
the Hilbert space to the infinite space. This point is espe- 
cially important for the soft modes, where amplitudes of 
vibrations are large and may exceed the range (maximal 
(a^)) of the physical space. 

We illustrate this problem in the Lipkin model where 
we know the physical space exactly. The HPT (|59p maps 
the angular momentum space {|JM)} onto the phonon 
space {\n)} (see Ref. [9]): 



\JM) 



\n = M + J), 



(70) 



where \n) is the eigenstate of A^ A. Since —J < M < J, 
we have <n<2J. By eq. ((511), 2J = n{n + 1) is just 
the valence particle number. 

Now we consider the possibility of diagonalizing eq. 
([62)1 in the infinite space. The negative A^^'^'> n^ / A term 
causes divergence. Thus we have two steps of approxi- 
mations: first, the A'^'^^V^/4 term can be neglected when 
diagonalizing eq. (15^ in the physical space {\n < 2J)}; 
second, the space can be increased to the infinite space 
{\n< -hex))}. 

The negative A'^''^V'*/4 term is smaller than the 77^/2 
term in the physical space (especially for the first few 
excited states), on the w^ > side of the critical point. 
Eqs. ((6T|) and ((Ml) give 



|A("4)| 



1 



'4J + 4kJ2 



1/k' - 8J' 



(71) 



The equality sign in eq. (|71|) holds at the critical point 
when K = Kf. = 1/(2J). On the a;^ > side 

{n\Tr^\n) _ (6n2 + 6n + 3)/(4w4) 6n^ + 6n + 5 

^ : ( ' ^) 



{2n+l)/{2v^) 



2n-|-l 



where the equality sign holds at the critical point when 



1/2. Consequently 



^A(04)^4/4^ 



72) 



2 8J 



(73) 



The upper limit of eq. ([75)1 is reached at the critical 
point for the state with the maximal number of phonons. 
We see that in the physical space the negative A(°^^7r''/4 
term does not reverse the order of states. For the first 
few excited states the upper limit in eq. ([75]) is actually 
much smaller, of the order 1/ J, because the upper limit 
in eq. (|72p is of the order 1. 

The space can be safely increased to the infinite space 
when J is large enough. The range of the physical 
space {n — 2J|a2|7^ — 2 J) ~ J increases linearly with 
J. The zero-point vibrations (a^) in the first few ex- 
cited states also increase, but much slower. On the 
a;2 > side, an upper limit is obtained when dropping 
the harmonic potential w^a^/2 in eq. ([5^ . in which case 
(q;2) ^ (A'^^"))~3 ~ Js. However, it is not justified when 



the collectivity is not so large, or if A^^") is numerically 
small (thus large zero-point vibrations, see Sec. IVB[) . 

We do a numerical example to illustrate the above two 
steps of approximations. The results for the first exci- 
tation energy Ei — Eq, at the critical point w^ = 0, are 
presented in Table HI In the last two lines eq. (|58)) is 
diagonalized directly in the {\JM)} space, where k takes 
the critical value corresponding to w^ = 0. In the last 
line the critical k is calculated by the RPA secular equa- 
tion ((66)) . with n = 1. In the second last line the critical 
K is calculated from 



1-4k2(J2 + J). 



(74) 



Eq. ([TUl is better than eq. (|551) because it is accurate 
not only in the leading order but also in the next order 
of 1/ J. 

The difference between line 1 and line 4 comes from 
neglecting higher orders in 1/J of A^™"^; between line 
1 and line 2 from neglecting the negative — gjx term; 
between line 2 and line 3 from increasing the space. We 
see that they agree quite well, and better for larger J. 
The difference between line 4 and line 5 is because the 
RPA secular equation is accurate in the leading order of 
1/J but not in the next order, which is the source of the 
biggest error in our method. 

In summary we argue that the existence of a finite 
physical boson space is general, in which the bosonic 
Hamiltonian should be diagonalized. This Hamiltonian 
may have "divergent-looking" terms [e.g. the negative 
^(04) ^gj-jjj ij^ gq_ ([g5])], which are indeed well-behaved in 
the finite physical space. 

However in general the exact physical space is un- 
known. Further approximations are needed if the mi- 
croscopically calculated (e.g. by GDM) bosonic Hamil- 
tonian is used to reproduce the spectrum of the origi- 
nal fermionic Hamiltonian. First, the "divergent" terms 
must be small and have little influence on the interested 
quantities, thus they can be dropped. Second, the inter- 
ested quantities must have reached convergence within 
the physical space, thus formally the bosonic Hamilto- 
nian (without the "divergent" terms) can be diagonalized 
in the infinite boson space. If the above two conditions 
are not satisfied, the bosonic Hamiltonian encounters se- 
rious difficulties or might be inapplicable in reproducing 
the correct spectrum. 



V. FACTORIZABLE FORCE MODEL 

Here we consider the factorizable force model where 
the GDM method provides approximate analytical re- 
sults. They will be compared with the exact results ob- 
tained by the shell model diagonalization. First we in- 
troduce a Hermitian multipole operator 



Q 



y^gi2a|a2. 



(75) 



For simplicity we assume qi2 is real; its hermiticity im- 
plies qi2 = (721- Furthermore, we assume that Q is time- 
even. The model Hamiltonian is 

H = ^ eia\ai + t X! (~'^9i4'?23 + K9i3924)aI4«3a4- (76) 

1 1234 

By definition of this model, the two-body part is different 
from 



-g-g = -:rXl(99)i2aI 



a2 



-- ^ (-Kgi4923 + K.q\zq2-C)a\a\a-ia^ ijl) 



1234 

by a one-body term. 



A. The GDM Method 

The mapping of Q is performed by substituting eq. 
HI) into eq. ([751): 



Q{R} = Tr{gi?} = 



Tr{gp} + Tr{gr^^"^}a + Tr{gr("^)}7r -f- Trjqr^^"^} 
+Tr{gr(«^)}^ + Tr{gr(")}i^ + 



Q 



(oo)_^Q(io)^^0-^ + Q(20); 



{a,7r} 



+g("2)^ + o. 2 



.,(78) 



where Q^"") = Trjgr^™")}. If n is odd, Q^"") vanish 
since we assume that Q is time-even. All Q^™") are real 
since Q is Hermitian. The self-consistent field becomes 



W^{i?}i2 = ^(-Kqi2'734 + K'?14'?32)^43 



34 



'«9l2 V" 934^43 = -m\lQ{R}i 



34 



(79) 



where we make the usual approximation keeping only 
the "coherent" summation. This is obvious in the har- 



„(io) 

43 



943; 



monic order, where the justification can be 
for higher orders this approximation is discussed in Ap- 
pendix [D] Substituting eq. (fTS]) into eq. (f79| we obtain 
the expansion of W{K\\2- 

Below we summarize the main results. Details includ- 
ing solutions for rj^™" are given in Appendix |E1 In the 
mean- field order we solve the HF equation ([T6l) : 



[/,p]l2=0, /l2=^12-Kg^°°^(Zl2. 



(80) 



Having in mind a spherical mean field, we assume that 
in the solution Q^''") = Tr{qp) = 0. Thus / and Z are 
the same, e\ = t\. 
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In the harmonic order the RPA secular equation ([50]) 
becomes: 



+2(.Q(i°))4 • { E E 



gi3[9, (t -q) ]32g2i 



i = -E 



\qi2\ "-12612 



t^ (ei2)2-a.2 
The normahzation condition (15^ becomes 



i = -(.g(-'fETfe\aS 



12 



[(ei2)^ 



(81) 



(82) 



ei#e2 63(7^62) 

E E E 

615^62 62/ ( = 62) 3 



612632 

'?12''?2'3932 «23921 
612(623)^ 



For higher orders we give the leading order expressions 
in oj^, following the procedure of Sec. IIIIl In the cubic 
order, eq. ([55|) becomes 



A(3") ^ (^Q(io))3 ( ^ 

Ely' 

+ EE 



[g. (f : g) ]i2g2i 

, 612 

ei#e2 

gi2g2i'9i'i '^12 



ei— e^/ 2 



(612)^ 



(83) 



where we have introduced notations for the weight factors 
(f : q)i2 = ^qi2- In eq. §^, the substitution of kQ^^o) 

by the leading order of eq. (|82|) gives Ac . Another 
equation ([54|) becomes 



A(i^) ^ -2 ^ +4(«g(i"))2A(30) ^ (J : ,),2,2i 



61^62 



[ 9, (f : 9) 112921 



(ei2)3 



+2(«Q(io))3 • ( 2 ^ 

ei7^e2 

E [g7(^ : g)]i2g2i _ y^ y^ gi2g2i'gi'i "-12 

(612)2 ^ 4" 

6151^62 ei=ei/ 2 



(612)^ 



.(84) 



_L^ V V [ g'(t ■ g) ]l2g21'gl'l 

3 ^ ^ (ei2)2 

ei=ei, e2(#ei) ^ '^^^ 

3 ^-^ f-/ . 612 

ei=ei/ 62(5^61) 

In eq. (|85|) there is an undetermined parameter A*^^^^. 



B. Two-Level Model 

Here the GDM method is compared with the ex- 
act diagonalization in a simple two-level model (see fig- 
ure4.jpg). The model has two s.p. levels with energies 
±i (the spacing is the energy unit), each with degener- 
acy n. There are N = n particles. The nonzero matrix 
elements of q are: vertical qi2 — q2i = 1, gii' = a/b for 
the nearest neighbors of the +2/ ~ 5 s.p. levels (the left- 
most and rightmost s.p. levels are also connected by a/6). 
Each s.p. level 1 is connected to only a few (three) other 
s.p. levels by qi2, thus the approximation in eq. (I79p is 
justified, as explained in Appendix |D] In summary the 
interaction has three parameters: overall strength k, and 
ratios a, b. 

In the mean-field order, Q^°^'> — X^i"-!?!! — since 
qii = 0. Hence s.p. energies e are the same as e. In 
the harmonic order, only the vertical qi2 = 921 = 1 ma- 
trix elements contribute. The RPA secular equation ([5T|) 
becomes 



Eq. dHU) determines Q^^'^-' as a function of A^^^^. 

Summarizing the results in this order: there are two 
undetermined parameters Ac and Qc ; Ai Ms fully 
determined; Qc and r^i2 1 ^ci2 ^'"^ determined as a 
function of Ac , r^-^^ is determined as a function of 
Ac and Qc .In the present model Qc , Qc play 
the role of the "undetermined overall factor" in w\2 \ 

,(02) 



1 - cj^ == 2Kn. 

Using eq. (|86p the normalization condition 
collective amplitude 



(Q(i°))2 = 2f7. 



(86) 
gives the 



(87) 



In the cubic anharmonicity, A 



(30) 



by eq. 



w{2 of Sec. IIIH respectively. 

The quartic potential term ()55p becomes 



there is no way to complete a three-body loop. In the 
quartic anharmonicity, Ac is calculated from eq. ([85^ : 



A(«)=2(^Q(i°))2(A(3"))2 ^(^: 5)^,5,^ 

_A(30)a(^2) 

+ (^.Q(io))3A(3°) . { 2 y ^^^^\-^\}^^^'^ 
[q, [^ ■■ q)]i2q2i , ^^ [q, (^ : g)]i2g2i 



Al^o) 



1 



[l-2(a-6)^ 



(88) 



E 

£15^62 



(612)= 



-+E 



ei7^e2 



-2EE 

e\—e-^f 2 



612 

»12gl2g21'gl'l 
(612)4 



The numerical diagonalization is done at $7 = 8, thus 
Kc = 1/16 and Ai*°' = [1 - 2(a - bf]/8. 

First we set the parameters a — b — Q. Fig. [T] shows 
the first excitation energy Ei — Eq as a function of k (for 
now ignore the two dot lines ^^a — b ~ 0.353" and "a ~ 
—b = 0.353"). As K increases to the critical value, the 
RPA frequency a; drops to zero, while Ei — Eq with the 
quartic potential term A*^4°) remains finite and agrees well 
with the exact results. This term restores the stability 
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of the system near the critical point. We emphasize that 

we have replaced A('*°) by A^''"^ = 1/8 in eq. gZ]), thus 
we are making a big mistake when uP' is large. However, 
it does not matter too much since in this region uP'o? jl 
dominates over A('*°)a'*/4. 

Next we consider the case of nonzero a and h. Both 
the exact E\ — Eq and our collective Hamiltonian are 
invariant under the change (a, 6) -^ {—a,—b), thus it is 
enough to consider only positive a. From the three lines 
of Fig. [T] "a = 6 = 0", ''a = b^ 0.353" and "a = -6 = 
0.353" , we see that the exact Ei — Eq depends on a — 6, 
but is almost independent of a + 6. This is in agreement 
with our collective Hamiltonian: lo is independent of a, 
b (in leading order of 1/fi); A^^^^ depends on a — 6 but 
not on a + 6. In the region of small k, the w^ potential 
term dominates thus Ei — Eq depends weakly on a — 6; 
whereas in the region of k ~ Kc, the A'-^'^-* potential term 
is important thus Ei — Eq depends relatively strongly on 
a — b. 

As a = — & increases, Ac decreases. At some point 
Ac becomes small numerically and ^ — I " ^ " no 



2 ' 4 ^ 

longer describes the behavior of the system near the crit- 
ical point. First, other anharmonic terms, suppressed 
by powers of l/fJ, may become important (see Appendix 
IF|) . Second, even if there are no other anharmonicities, 
the description breaks down because the increasing zero- 
point- vibrations (a^) will exceed the range of the physi- 
cal space, as discussed in Sec. IIVDI The current model 
has a larger vibrational amplitude than the Lipkin model 
due to a smaller Ac (^ 1/Vt verses 4/J7). Fig. [5] shows 
El — Eq as a, function of the parameter a = — 6 at the 
critical point to — [k, = Kc = 1/16). Ei — Eq depends 



on the space in which we diagonalize ^ 



A(40)^ 



Unlike 



2 ' 4 ^ 

in the Lipkin model, we do not know a priori what the 
physical space is in the current model. But it should be 
similar to that of the Lipkin model with 8 particles. Thus 
we choose rtmax — 8 for both two finite spaces, each with 
a reasonable u of eq. (|60l) . When a = — 6 is small, say, 
less than 0.2, Ei — Eq of different spaces are close and 
all follow the trend of the exact Ei — Eq. When a = —b 
is large, Ei — Eq of different spaces differ substantially, 
implying that (a^) has reached the edge of the physical 
space, thus the bosonic approach becomes invalid. If in 
the current model we increase the collectivity 17 = A^, it is 
expected that Ei — Eq from the GDM method will agree 
with the exact Ei — Eq up to a larger value of a = — &. 

In summary, near the critical point w^ ~ 0, the next 
even potential term A^""^' dominates the dynamics of the 
system, provided it is positive and large. A*-'*^"-' should 
be large enough such that other anharmonicities were 
negligible, and zero-point vibrations (a^) were within the 
finite physical boson space. A larger collectivity factor 
51 helps both, since other anharmonicities are suppressed 
by powers oi fl^^ (see Appendix IF|) . and the range of the 
physical space grows as f2. 



VI. REALISTIC NUCLEAR APPLICATION 

There are three complications in realistic applications 
of the GDM method. A realistic nucleus has two kinds of 
fermions; symmetries, e.g. rotational invariance, need to 
be respected; pairing correlations should be considered. 

As in the BCS theory we substitute the original system 
by a grand-canonical ensemble, in which the chemical 
potential is fixed by the average particle number of the 
ground state in the mean-field order. In this case we need 
to consider e.o.m. of not only 0201 but also 0201. A good 
treatment of the superfluid ground state, on top of which 
collective excitations are formed, is essential. 

The collective mode operators a\^ , tt^^ have quantum 
numbers corresponding to symmetries of the Hamilto- 
nian. In this section we keep only the quadrupole mode 
which is the most important one at low energy. The case 
of interacting modes (quadrupole and octupole) is dis- 
cussed briefiy in Appendix [Lj 

This section is a straightforward generalization of Sec. 
im The details of the derivation are given in Appendix 



A. Preparation 



The microscopic fermionic Hamiltonian for the canon- 
ical ensemble is still given by eq. ([T|): we include the 
—IJ,N term in Z, and the s.p. index 1, 2 ... can run over 
protons and neutrons. Isospin may not be conserved for 
some effective interactions. We do not write V in the 
form ^Jij2),(j3j4); ^12 and ^1234 carry all the symmetries 
of H implicitly. 

Now the reference state |$) does not have definite par- 
ticle number, 

($|aja2|$) = P21, ($|aia2|$) = K21. (89) 

K is the pair correlator [28[ . Also we need two generalized 
density matrix operators 



,t. 



(90) 



i?i2 = a^ai, K12 = a2ai, 
and two self-consistent field operators 

Ty{i?}l2 = ^yi432i?34, f{R} = Z + W{R}, (91 



34 



A{if}i2 = |^yi234if43. (92) 



34 
T fT 



It will be convenient to introduce {R , / are transpose) 



D{R,K}^[ ^, j^^^ 



S{R,K\= (^At{Al -r{R}) 



(93) 



The collective mode operators a^ , ttI carry quantum 
numbers of angular momentum A, its projection /i, and 
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parity (— )'^. The coordinate a\ is time-even, and the 

momentum ttI is time-odd. Their Hermitian properties 
are 

"U - (-)'"''«A_^, 4, = (-)'"^^A-p. (94) 
The commutation relation is given by 

["L'^AV'] =« ^AA"5^M'- (95) 

Here we consider only the quadrupole mode A = 2, and 
drop the label A. 

The collective Hamiltonian replacing eq. ^ should be 
written with correct vector coupling of the operators: 

H^Eo + Y\/5(a X af„ + -\/5(7r x nf„ 
+ —--V5{{a X a)2, a}[] + ^— \/5{a, (tt x t:)^}'^ 



A(40) 

H — — \/5((a X af x (a x a)") 

A(o4) 

H — V5((7r X 7r)° x (vr x 7r)°) 

(22) 







+ E ^V5{(axa)^(^x^)^}«.(96) 

L=0,2,4 

7f is Hermitian, time-even, invariant under rotation and 
inversion. 



+ [5(10)^^(01)]L„[^(10)^^(01)]L_(^02) 

In the quartic order: 

-|-' E <.'-7.^,.-3^(-)^A(-^<^-7.\ 

1=0,2,4 

+1 E i[sr\D^'x~[Dr,s''X)-^t.Am 



{(qT X qT)', ^T}i./4. ; = 0, 2, 4; L = 0, 1, 2, 3, 4, 5, 6 



2i 



-2^(-)-.;^ E D^Z-9t' + ^S,,.D 



(30) 
Lfj. 



+2zA(i^'<«'.g,^,+.5^.Ari3r' 

^j5(20)^^(01)jL_j^(2a)^^(01)jL 

-2 E ([5r',^''°']^-[-Dr'>^'"']^)-5u', (104) 

('=0,1,2,3,4 



B. Equations of Motion in the Collective Band 

Following the same procedure as in Sec. [Hi we find 
e.o.m. replacing those in Sec. IH BIIH El Matrices d(™"), 
Simn) j^j.g coefficients of expanding D{R,K}, S{R,K} 
over collective operators a', ttJ^. In the mean- field order 
we obtain the HFB equation 

[5{p,4,Z?{p,4] = [5("o),Z?("o)]=0. (97) 

In the harmonic order we obtain the QRPA equations 

4: zi^(i°) = [5(™),i?f)] + [5(°i),i?(°")], (98) 

al : - tu'Dl^'^ = [S'^'"\dI!"^] + [S(^"\d^°% (99) 



In the cubic order: 

(atxat);;/2,X^ 0,2,4: 

-2tLu^D[''J ~ 2*<5L,2A(3")i?(^oi) = 
[^(oo)^^(20)j ^ [5(^20)^^(00)] ^2 [5(i"),Z?(i")]^, (100) 
(7rtx^t)i^/2,X^ 0,2,4: 

[5(00)^ ^(02)j ^ [5(02)^^(00)] ^2 [5("i),Z?("i)]^, (101) 



+2z(-)^ E i?g^^y^-i5,,A(-)i3 



K, (tt^ x ■K^YJ^/i; l = 0,2A;L = 0, 1, 2, 3, 4, 5, 6 : 

2*c / x-L 2r-i(03) 4i 2r-,(03). \L-lr L 

(i2)n(ii) 

Ljj. 

+2zH^A(^^'<^ 5.", - ^ 5,. Ar i^r 

^j^(10)^^(02)jL__j^(10)^^(02)]^ 

+2 E ([4"'.^^"']m - [D'-^'\s''''t){~)'-'9t', (105) 

;'=0,1,2,3,4 



{(7rtx^tyL^^t|L^e^^^0^2,3,4,6: 

— ( ^^ "S^ d(^^) -v-^ ^a(^2'd(°2) -v-^ 

;=0,2,4 

+3J5i2A(o^'75(/») . 7o^r,' = [S(»«\i3r'] + [Sr',Z?("»'] 

1=0,2,4 

The numerical coefficients 7^^;, and 5;^;, are defined by 
{(a X «)',«};: =7i''r ■{(« X a)'',a}^^, (7,^, = 1),(107) 

;=0,2,4 
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Values of -f^i, and 5/";, are given in Appendix ID II 
in {{a^ X a^)'^,a^}f'/6 is the choice of basis, different 
choices of II do not influence results. 

There exists a relation involving cubic and quartic an- 
harmonicities, replacing eq. (|331). Setting L = 2, keep- 
ing only (i5/5)(30/2i/i2/03) and A(40), Ap 2)^ a(04) terms, 
/;^x eq. ([Ml) +5W^Ei=o,2.4/;X eq- dTOH) gives 

i=0,2,4 

= [ ^<°°^ ,(/..• /?f ' + i-^ E /' • ^S^ ) ] 

;=0,2,4 

1=0,2,4: 

+3,A(40)^(01) . f^ + ^^2 ^ ^^ . A(22)^(01) ^ (^og) 

i=0,2,4 

lJ2i=o,2,4fi'x eq- dini +wVi^x eq. dM]) gives 

;=0,2,4 
i=0,2,4 



where f\{r) is real. The definition of eq. (J112I) differs 
from the "usual" one in two aspects: a factor i^ is in- 
cluded, and ql ^ Y\^ instead of q\fj_, thus ql creates 
projection fi. The Hermitian properties are 

9L = H'^^gAp, Q{, = {-)'-''Qx,. (113) 

The pairing operators P and P^ are defined by 



pt^ly-alal P=i^ 



2 / ."-l^-T' 

1 



2/l^«i«i' 



(114) 



(^ E /.-^S^-v^.-^ 



(03) \ ^(00) 



/=0,2,4 



-- Y^ fi- A\^^'>D'i}°^ -3iuj^A'^^'^^Dl^°^ ■ fo + ...{no) neutron, and proton-neutron (Kp 



;=0,2,4 



where fi is defined in eq. (jlTj). A solvability condition 
exists because the variable parts of eqs. (|109p and (jllOp 
have the same structure as the QRPA equations ((98)) and 



where 1 is the time-reversed s.p. level of 1. P has angular 
momentum and positive parity. P -I- P''^ is time-even, 
P — P^ is time-odd. 

The quadrupole plus pairing Hamiltonian is (dropping 
A = 2) 

iJ = ^(ei - n)a\ai " j E ^WiO'~2"'2 
1 12 

+ 4 E E(~'* <?ii4<?A'23 + li qli3q^24:)a\ala3a4. (115) 

1234 M 

Approximately, this Hamiltonian can be written as iJ ~ 
J2i (ei - M)alai - GP^P - \k ^^ Q\Q^. The difference 
is in a one-body term originating from the Q ■ Q part. H 
is Hermitian and time-even, which implies real G, k, ei = 
ej. In a realistic nucleus there are protons and neutrons; 
formally we can still use eq. (|115p if the quadrupole 
force strengths are the same for proton-proton, neutron- 

n = Kp„ = k), while 
remembering the pairing is treated for protons and neu- 
trons separately {Gp ^ Gn). We will assume this is the 
case. 



Following the procedure in Sec. IHIl we can obtain 
expressions of Ac and Ac .In the next section we do 
this explicitly for the quadrupole plus pairing model. 



VII. QUADRUPOLE PLUS PAIRING MODEL 

In this section the GDM method is applied to the 
quadrup ole p lus pairing Hamiltonian. As was understood 
long ago [23, [sol, this model combines the most important 
nuclear collective phenomena in particle-particle (pair- 
ing) and particle-hole (quadrupole mode) channels. The 
approximate analytical results of the GDM method are 
compared below with the exact results of the shell model 
diagonalization. The operator of multiple moment is de- 
fined as 



A. The GDM Method 

1. BCS 

In the pairing plus quadrupole model the HFB equa- 
tion (|97)) becomes the BCS equation: 



1 ^ 

ei ^ ei - 11 - G{vif , (117) 

E^ = v/(ei)2 + (A)2, (118) 

(«i)'-^(l + ^), (^i)' = ^(l-|^)' (119) 

N^yiv^f. (120) 



1 



QUR} - Tr{ql^R} = J2 9Li2«l«2, 



(111) 



qi^ = h{r)-t^Y^^{9,cb), (112) 



BCS amplitudes ui = ui, vi = vi are real. Pairing 
energy A is a real number, not to be confused with the 
field AjK} in eq. (|92p that is an operator matrix. Ei 
is the quasiparticle energy. The chemical potential fi is 
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fixed by eq. (|120p . The gap equation (|116|) lias a non- 
trivial solution A > only if G is greater than its critical 
value Gc |30|- For convenience we introduce: 

(2) _ {U1V2+U2V1) _ 

V12 = (P, _L P„\2 g^l2, ry^^lZ = («1W2 - 1)11^2)9^12 -(121) 



iEi + E2y 



2. QRPA 



The QRPA secular equation corresponding to eq. ([5T|) 
is given by 



1 = K 



12 



(£'1 + £'2) \S,fj,2l\ 

{Ei+E2Y-0J^' 



(122) 



The solution w^ is independent of /i. Results in the form 
of reduced matrix elements are given in Appendix|Kl The 
normalization condition corresponding to eq. (|82p is 



1=(«0(1"))2 ^,(^l+^2)|e,.21P 



(123) 



z-.[i^E,+E2Y-w^r 

3. Cubic Anharmonicity 
The cubic anharmonicity corresponding to eq. 
A(30)=3(Kg(i"))3 J2 

(2ji + l)(2j2 + l)(2j3 + l) f 2 2 2 
5 \ ii J2 33 

. <^(i)t ^(i)t t Q241 

^||12 ^||23 ''1131' V-l-^"iJ 

where O,^^ = {niji\\^^^'^'^\\n2J2) is the reduced matrix 
element, the convention for which is given in Appendix 
|j] rii combines all other quantum numbers specifying a 
s.p. level, except j'l. 

We give the expression of p(^") which will appear in 
a(^°): 



F(2").[l-G^(2ji + 1) 



[(«i)'-(«i) 



212 



riiji 



4£i 



.Q(io))2 E (2j-i + 1)-V27^TT{| I 

(Ui)2 - {vif J (i)t 



nijin2j2 



2 2 
h 



9„ „ ;^(l)t;,(l)t 
2-UlWl -41112 41121 



£1 



•'?ii2eii2;' -(125) 



p(20) jg divergent when G is greater than but close to Gc- 
In this region of the pairing phase transition, A is small, 
and F(20) _ 1/^ The GDM + BCS method is not valid 
in this region: in the mean-field order the BCS solution 
already fails, as is well known. 



4- Quartic Anharmonicity 
The quartic anharmonicity corresponding to eq. ([85 



IS 



/o . A(*°) = -2 /2 ■ («Q(i°')^(A(^°')^ Tr 






(5)t 



-Tr 



H^'''}r'e'' +2Tr 






+Tr 



{v,eY,^'e'' 



-/o • G p(^°)(.Q(^'")^ ^(^.1., + U2.2) e^i^e^i^ 

12 

+/o-GP(^")(A.Q(i'")^ 
^(ui)2-(i;i)^ , (^2)^-(t.2)%,(i) t 



■Ei^ 



+ 



"Ivi2''m21 



2Si 2£'2 

+ (.Q("'))^ E .A-(Tr[u,(?«xC<^')'}r'e<^^'' 

;=0,2,4 ^ 



-Tr 



{^(i),{^,^(i)y,(i)}^^2^t 



(126) 



where {?7,C^^^}ii^^ = {'7,?'^^}'i2/(^i + -£^2). There is an 
undetermined parameter A^^^^ in eq. p26p . Values of 
numerical factors /; are given in Appendix HI 

B. Comparison with Exact Results 

We compare the results of our method in a semi- 
realistic model with those of NuShellX [3l|. There are 
10 fermions of one kind and four s.p. levels with ener- 
gies: 



s.p. levels 
e (MeV) 



Ipi 
-0.1 



0/1 

o.o" 



IpI 

1.0 



0/i 

1.1 



We take the radial wavefunctions to be harmonic os- 
cillator ones. In eq. (|112p we take /(r) to be r^ so 
qjj^ = —•r^Y2fj,{9,(j)). For convenience we make qj^ dimen- 
sionless by combining its original dimension with k (see 
the end of Appendix [j|. The model space is similar to the 
realistic p/-shell, but the Ip^ and lp| levels are inverted 
to increase collectivity: in the current case the q matrix 
elements ((7i„i q/s and q^^ a g^i) are large between the 
s.p. levels above and below the Fermi surface. 

We did a set of calculations with increasing pairing 
strength G. At each value of G, the strength k of the 
Q ■ Q force is taken to be at the critical value Kc such that 
the RPA frequency w^ = 0. The results are summarized 
in Table |TT1 For clarity, we draw the last three lines of 
Table M as Fig. |31 The coefficient At^o) in Table M 
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is calculated by eq. (|126p setting A^^^^ — (dropping 
the — /2 • A'^^^A'-^^) term). A non-zero term A^^^^ in its 
reasonable range does not influence A^^*^' much, since in 
the current model A^^^'' is small due to the approximate 
symmetry with respect to the Fermi surface (see Table 
HI)) . Then "GDM i?2+" is calculated by diagonahzing eq. 
dSH), setting A(i2) ^ A^O'*) = A^^^^ = (w^ = since k 
takes its critical value). 

The critical value of the pairing strength Gc is around 
0.11 - 0.12 MeV. When G < Gc, the BCS solution A = 
0, and n can be anywhere between Eq^ ? = and e^^ 3 — 

1.0 MeV. We checked that in this case our results (I122f 
I126P do not depend on the choice of /i. In Table lU we fix 
fi at 0.5 MeV. In the region where G is greater than but 
close to Gc, our method is invalid as discussed under eq. 
(IT251) . This is iUustrated in Fig. |3]by the 'kink' on the 
"GDM E2+" curve near G - 0.12. 

In Fig. [3] "exact E2+" and "exact £^4+" are the exact 
results by NuShellX. At G = the first excited state is 
4+ instead of 2+. In this case the 4+ state is a single- 
particle excitation from 0/| to 0/|; the 2+ state is a 
collective state with approximately half holes in Ip-^ and 
0/| levels, half particles in lp| and 0/|. As G increases, 
the collective 2+ state becomes the first excited state. 
When G is large enough, A dominates over the original 
s.p. spacing e, and the results become stable. As an 
example, at G = 0.30 MeV, the quasiparticle continuum 
starts at ~ 3.5 MeV; from the second excited state 0+ at 
3.506 MeV to 4.153 MeV there are 15 states with J^ = 
0+,2+,4+,6+. The first excited state 2+ at 2.438 MeV 
should be identified as a collective state, stabilized at 
around sixty percents within the gap. 

It is seen in Fig. [3] that "GDM E2+" agrees well with 
the exact result "exact E2+" in general. On the G < Gc 
side, our E2+ does increase with G although not rapidly 
enough. On the G > Gc side, when A is not too small, 
the agreement is very good. 

In summary, this section shows the potential of the 
GDM method in doing realistic calculations. In medium 
and heavy nuclei the pairing gap ~ 2A ~ 2 MeV, the 
critical region is approximately bounded by |a;| < 1 MeV. 
Nuclei on the uj^ < side are gamma-unstable. On the 
u!^ > side, the whole region can be calculated as in Fig. 
[T] [explained in the paragraph under eq. 



VIII. CONCLUSIONS 

The GDM method is promising in solving the long- 
standing problem: constructing the collective bosonic 
Hamiltonian microscopically. The procedure is straight- 
forward and consistent. Results of the lowest orders, 
the well-known HFB and QRPA equations, give us con- 
fidence to proceed to higher order anharmonicities. The 
anharmonicities are important as the harmonic potential 
Lo^c? jl becomes small or negative when going away from 



closed shells. The GDM method provides a unified de- 
scription of different collective phenomena, including soft 
vibrational modes of large amplitudes, gamma-unstable 
potential and transition to static deformation. It maps 
the exact fermionic e.o.m. onto the dynamics generated 
by approximate collective operators. Here we used the 
phonon-like operators; other possibilities include rota- 
tional dynamics and the dynamics corresponding to the 
symplectic symmetry or other group-theoretical models. 
In such cases the GDM expansion should be based on the 
group generators. 

Sec. |TT] discusses the general procedure of the GDM 
method. In each order, a set of coupled linear equations 
is solved in terms of lower order results. At the current 
stage the anharmonicities are not completely fixed; we 
find only one relation (I44p involving the cubic and quar- 
tic anharmonicities, appearing in the third order as a 
solvability condition. In Sec. Illll it is shown that around 
the critical point w^ ss 0, we are able to determine the 
cubic potential A^'^'^^ (|55|) and the quartic potential A^'"') 
([55| . A*^^") should be responsible for restoring the stabil- 
ity of the system near the critical point, if it is positive 
and large. This idea is then tested in three models of in- 
creasing complexity: the Lipkin model (Sec. lIVp . model 
with factorizable forces (Sec. |V|, and the quadrupole 
plus pairing model (Sec. IVIip . The GDM method is only 
responsible for calculating A^™"'; other conditions are 
needed if the resultant bosonic Hamiltonian is used to 
reproduce the spectrum of the original fermionic Hamil- 
tonian, as discussed in the last two paragraphs of Sec. 
IIVDI If these conditions are not fulfilled, the approach 
of effective bosonic Hamiltonian encounters serious diffi- 
culties. The conditions for quartic potential (A'^'*°^a''/4) 
dominance near the critical point are discussed in the last 
paragraph of Sec. IV Bl 

Calculations for realistic nuclei are in progress. How- 
ever, the pairing correlations need to be treated better 
than in the BGS framework, because anharmonicities 
are sensitive to the occupation numbers (u\, v\) of the 
superfluid ground state. Unlike the QRPA secular equa- 
tion (jl22p . where terms in the summation contribute 
coherently, in the expressions of anharmonicities (|124|) 
and (|126p different terms may cancel. A(^°) and A^ 
depend on the balancing above and below the Fermi 
surface, thus they are sensitive to the occupation num- 
bers (mi, v\). Work is also in progress about the role of 
A(30) on deformation, as well as the quadrupole-octupole 
coupling in the presence of a low-lying octupole mode. 
The realistic effective interactions (better than the 
quadrupole plus pairing Hamiltonian) are to be used in 
the calculation. The present paper sets the scene for the 
GDM method in the sense that it is seen explicitly there 
are no contradictions in the solutions (Sec. HIl and IVTI . 
although at the current stage we find only one constraint 
(|1^ on the anharmonicities. New constraints, if found, 
would fix the anharmonicities completely. 
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Appendix A: Three-Body Force 

It is straightforward to include three-body forces in the 
formulation. The microscopic Hamiltonian ([1} includes 
a new (anti-symmetrized) term 

H^^'> = ^ X! G'i23456aI4"3fl4a5a6- (Al) 

123456 

Under the definition 



G{i?}l256 = 2^ Gi23456^43 
34 

G{R, R}i4 = 2_^ G{i?}i234i?32 



(A2) 
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the normal ordering Hamiltonian (jlOp acquires new 
terms, 

(<l>|ff(3)|<i>)^iTr[G{p,p}p], 
/^'^-^G{p,p}, F(3)=GM, (A3) 

and a term ;^ X]i23456 ^i23456A^[aia2a3a4a5a6]. In the 



e.o.m. (J12p / and V are replaced by the new ones in- 
cluding f^^^ and V''^^ {W{R} is calculated from the new 
V), and there are two additional terms: 

+ -[G{R'^^},p]i2 + T^ 51 iGi34567N[alalala5aea7] 



34567 



G765432A^[a7aJa5a4a3ai]),(A4) 



where 



we have the following identities in the full space: 

[i?i2,a] = [.T,i?]i2, (B3) 

[Ri2,tt]^[p,R]i2, (B4) 

where [a;, i?]i2 = J^ai^is^s'^ ~ Ri3X32)- Similarly to the 
manipulation of eq. p^ . we project eqs. (JB3IB4I) onto 
the collective subspace. Since a and it are collective op- 
erators, we can substitute R by its boson expansion (fT4l) . 
After calculating commutators on the l.h.s. , we equate 

2 

coefficients of the same phonon structure: 1, a, n, ^ ... 
Eq. (US]) gives 

-ir\^2^' ^[x,p]i2 = Xi2{n2~ni), (B5) 

(B6) 



-ir(02) = b,r(0i)l 



ir 



(11) 



x,r 



(10)1 



Eq. (|B4|) gives 



^(20) 



"12°^ = [P,P]l2 



= Pi2("2 -ni), 
•(11) = [p.r(°^)] 



(B7) 
(B8) 



Only the rii ^ 712 matrix elements of x and p are de- 
termined from eqs. (|B5[) and (|B7I) . Higher order ex- 
pressions (jB6p and (jBSp are approximate, saying that 
7.(20/11/02) g^j,g completely fixed by the the harmonic or- 
der solutions. In fact the two expressions of r^n) are not 
consistent with each other. These defects are due to the 
neglected many-body components in eqs. (|B1|) and (JB2|, 
as explained in Appendix [C] The approximate expres- 
sions (jB6|) and (|B8|) are used below to derive expressions 



01 r^234 m terms 01 r j^j 

In the full space we also have 

[A^[a4a3a2ai],a] = 

- 2J X25 N[alalaia5] + 2_] x 15 N[alala2a5] 

5 5 

- X] a;64A^[a64'^2ai] + ^ xe3N[alala2ai] 

6 6 

-hir^4^^7V[4ai] - ir^°^^^ N[alai] 
-tr^H'^ N[ala2] + tr^li'^ N [aia2l (B9) 



G{R }i6 = 2_^ Gi23456-^[a2"3'^4a5]- (A5) where we have used eq. (jBSp . Similarly 
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Formally the HF and RPA equations are the same as 
before, replacing / and W{i?} by the new ones. 



Appendix B: Saturation Principle for Section ITll 

Keeping only one-body terms in a and tt, 

a^^xi2R2i=^xi2a\a2, (Bl) 

12 12 

TT = ^Pl2i?21 = ^Pi2ala2, (B2) 



[N[alala2ai],n] = 

'^P25N[alalaia5] + y^Pi5A^[a4"3Q2a5] 

5 5 

■ ^P64N[alala2ai] +^pe3N[alala2ai] 

6 6 

-ir^^f^ N[alai] + ir^^f^ N[alai] 
+ir[\°^ N[ala2] - ir[f^ N[ala2], (BIO) 



12 



12 



where we have used eq. (jB7p . Again we project eqs. (|B9[) 
and (jBlOp onto the collective subspace, then substitute 
the expansions (fTi|) and (jlSp . Both the l.h.s. and the 
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r.h.s. have no constant terms. This justifies the assump- 
tion under eq. (1151) : terms hnear in a and tt are absent 
in the expansion (1151) oi N[a\a^a2ai], since they generate 
constant terms in the Lh.s. of eqs. (IB9[) and (JBIOI) . The 
a and n terms of eqs. (IB9|) and (|B10|) give 



r 



(10)^(10) ^(10)^(10) ^(10)^(10) .^(10)^(10) 



14 ' 23 



13 '24 



'24 '13 



'23 '14 



^(20) 
'1234 

(02) _ (01) (01) _ (01) (01) _ (01) (01) (01) (01) .„,„) 
'1234 — '14 '23 'l3 '24 '24 ' 13 +'23 ' 14 A^^^) 
„(11) _^(10)„(01) _„(10) (01) _ (10) (01) (10) (01) .„,oN 
'1234 — ' 14 '23 '13 '24 '24 '13 ~r'23 '14 -1,010) 

We mention that eq. (|B9I) and (|B10p give the same ex- 



pression of 



.(11) 

1234 



(|BT3|) . Using eqs. ^, ^ and ([BTH 



lBT3l) the a^/2, {a,n}/2, 7rV2 terms of eqs. (|B9| and 
(BTol) give 



J30) 



;('- 



(20) (10) 
14 '23 



p(20)^{10) 
13 '24 



-,{20)„(10) 
'24 '13 



+ r. 



(20) (10) 
23 ' 14 



.{03) _ 3, (02)^(01) ^(02)^(01) ^(02)^(01) .^(02)^(01) 



.(21) 



+ rL, 'r] 



),(B14) 
),(B15) 



J20) (01) 
'14 '23 



.(20)^(01) 



.(20)^(01) 



'24 '13 



^(20)„(01) 
'23 '14 



,2^(11)^(10) _2r'ii)r'i°' 

T^^'l4 '23 ^'13 '24 



2r, 



(11)^(10) 



+ 2r. 



(11)^(10) 
23 '14 1 



(B16) 



„(12) _ „(10) (02) _ (10) (02) _ (10) (02) (10) (02) 

'1234 ^'14 '23 '13 '24 '24 '13 "r'23 '14 

,9^(11)^(01) _ 2r<liV'°i> - 2r<llV(''l' 4- 2r<llV(''l' (Rl?) 
+^'l4 '23 ^'13 '24 ^'24 '13 ^^^'23 '14 -VlJl-'J 



Eq. (|B9p and eq. (IB10[) give the same expression of rl 



(21) 
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(IBT6)) and rl^H (|BT7)) . The resuhs (|BimB17p generahze 
the so-cahed hnearization of e.o.m. method, 

0^030201 -^ Piiala2 - Pi3ala2 - P2ialai + p23a4ai-(B18) 

The normahzation of the RPA sohition A^^\ r'^i) is 
determined by the commutator [a, n] = i. Under the 
one-body assumption (JB1[) and (jB2[) . 



i = [a,7r] — Tt{[x,p]R} 
= Tt{[x,p]p} + Tr{[a;,p]r(i°)}a + ... (B19) 



The constant term of eq. (jB19p gives eq. ([5^ . The 
higher order terms of eq. (|B19p should vanish, as dis- 
cussed in Appendix [C] 



Appendix C: Many-Body Mode Operators 

Outside the harmonic regime the mode operators a 
and TT have many-body components. Here we write down 
the results for a only, tt is treated similarly. The structure 
of a replacing eq. (IBip is 



a = 2J 2;i2a{a2 + t ^ xi23iN[a\ala3a4] 

12 1234 

-t-- ^ xi23456N[alalala4a5ae] + ■ . ■ , 



(CI) 



123456 



where X1234 and X123456 are anti-symmetrized structure 
coefficients. The saturation principle replacing eq. (jBSp 



[N[alai],a] = [x,p]i2 + [x,R%2 + [i{R^},p]i2 
+ - ^(:Ei345A^[44«4a5] - X5432N[alala3ai]) 



345 
,(B11) +[i{i?2'^},p]i2 + i ^(i 



il34567A^[a2«3°^4"5a6a7 



34567 



i765432A^[a7a605a4a3ai]) -l- . . . ,(C2) 



where 



23 



~S:{R''^} 



i{i?^}l4 = y^Xi234A^[4"3] 

Cl 
2345 
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22 il23456A^[a2«3"4a5] 



(C3) 
(C4) 



Comparing coefficients of the same phonon structure we 
obtain 



zr 



(01) 



[x,p], 



and 



and 



,(02) 



.("1)1 



ir 



-ir 



(21) 

12 



(11) 



\x,r 



(10)1 



[x("i),p], 
[i^''\pl 



(C5) 

(C6) 

(C7) 



b,r(20)]i2 + 2f5;(i"),r(io)l 



+ [x(20)+i^^"\p]l2, 


(C8) 


-zr(i/) = [x,r(")]i2 + [i(^°\r(0i)]i2 




+ [x(0i),r(i°)]i2 + [:£(") +^^"\p]i2, 


(C9) 


2zr(f = [x,r(02)]i2-f2[5;(oi),r(oi)]i2 




+ [^(-)+5^''^\p],2. 


(CIO) 



From eqs. (JCSIICIOP the structure coefficients x, x and x 
of a are determined by the e.o.m. solutions r(™") , order 
by order. For self-consistency, substituting them into eq. 
(|C1|) should give a, 



a = Trjxp} + Tr{a;r(i°)}a + Tr{ 



.(01) 



}^ 



-(Tr{a;r(2")} + Tr{x(i°)r(i°)})^ 



-(Tr{a;r("^^} + Tr{x("iV("i^}) 



(Cll) 



which means that all other coefficients vanish, except 
Tr{a;r(i°)} = 1. Trjxp} — implies that diagonal matrix 
elements xn = 0. TrjxrC^i)} = is satisfied identically 
by eq. (jCSp . Tr{xr(ii*)} = 1 is identical to the normal- 
ization condition ([5^ . For higher order coefficients in eq. 
(jClip . some are identically zero, e.g. the 7r^/2 coefficient 
by eqs. (|C5p and (|C6p : some impose new constraints, e.g. 
the vanishing of the a^ /2 coefficient implies 



.(11) 



Tr{a;r(-'i^} -f Tr{r(iiV} - iTv{[x, r^'''>]p) = (C12) 



(10)1 
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In the Lipkin model we have checked that these con- 
straints are satisfied identically, up to the a'^ , {a^,7r}, 
{a,7r^} and tt'^ terms. 

These many-body components should be kept in mind 
if we want to compare the bosonic wavefunction with the 
shell-model wavefunction. 



The leading order of eq. (IE1[) is 



The leading order of the RPA secular equation (ISTj) is 



Appendix D: Coherent Summation 

The factorizable force model has an analytical solution 
only if we neglect the "incoherent" terms in eq. ([7^ . as 
is usually assumed in such models. Here we consider 
its justification beyond the harmonic order. The exact 
expression of W{R}i2 is 



[van] 



34 



(-^512534 + K,quq:i2)r. 



X^ 



Kgi2 2^ 934^3 
34 



(mn) 



43 

Einin] 
34 



An observable is given by a trace of w^™" with some 
operator(s) t: 

O - Triiu-^"")] = 
^K ■ Tr[ig] • Tr[gr("")] + k ■ Tr[gigr(™")]. (Dl) 

Quite generally, operator qi2 has the following property: 
for a given s.p. level 1, gi2 essentially vanishes except for 
a few s.p. level 2. For the realistic quadrupole moment 
operator q^ = r'^Yi^i, it is ensured by the selection rules 
with respect to r^, L = 2 and fj,. If gi2 has the above 
property, a trace grows linearly with the collectivity fac- 
tor rj, independently of the number of operators q inside. 
Hence in eq. (|D1I) the incoherent sum is smaller by a fac- 
tor of 1/57 than the coherent one. The approximation of 
keeping only coherent terms is valid when the collectivity 
ri is large. 



Appendix E: Details of Factorizable Force Model 



Here we supply the details for Sec. IV Al In the har- 
monic order we solve the RPA equation. The formal 
solutions ([^ and ([^ become 



Jio) _ 



-kQ(io) 



912 



(01) _ 



(ei2)2 -uj' 



j"-i2ei2, r{^ = 



-kQ(io) 



912 . 



(ei2)2-w^ 



mi2.(El) 



Then Q("i) = Tr{qr^^'^'>} = 0, as it should be. From 
Q(io) ^ Trjgrfi"'} 7^ we obtain the RPA secular equa- 
tion (j8ip . The ni 7^ 71,2 matrix elements of x and p are 
given by eqs. (IBSP and (|B7|) : 



a;i2 



(ei2)2-^2 



P12 



kQ(io) 



912 



(ei2)^ 



iei2, (ni 7^ n2).(E2) 



i = -E 



Igi2r»'i2 

ei2 



(E4) 



The leading order of the normalization condition (|5^ is 

|2, 



1 = - 



(,Q(10))2^kl^. 



(E5) 



12 



In the cubic order, the ei = ei' matrix elements are 
given by eq. ([55)) : 



(f ^-2(.Q(i°)f^ 



(10)n2 Y^ gl2g21' '•^12 



(ei2)^ 



^(02) 
'll' 



V (ei2)* 



'll' "■ 



(E6) 

(E7) 
(E8) 



The ei 7^ 62 matrix elements are determined from eqs. 
(13311351) . 

-2zA(^o)r(f =ei2r[f +«Q(2")ni29i2 

-2«Q(io)[9,r(io)]i2, (E9) 

irif + .A(i2),(io) ^ ^^^^(11) _ «,Q(io)[,,,(oi)]^^^(Eio) 

2zr(i/) - zA(i2),(oi) ^ ^^^^(02) ^ ^Q(02)„^^^^^^ (Ell) 

with the solution (ei 7^ 62) 

.(f = -.Q(20)(!^ : ,)i2 - 2.Q(i°)A(3°)(4 : q),, 
e e^ 

_2(«;Q(io))2kiil^,(E12) 



ei2 
n 



.(i/)^_,«Q(20)(4:,), 



z..g(i°)[2A(30)(4:9)i,+A(i2)(4:g) 



-z(«:Q(io))2{2 



[9,(?:<z)]i2 , [g,(^:g)] 



12 



„(02) 



(ei2) 
,n 

' e 



ei2 



},(E13) 



^Q''''\-:q)l2 + 2nQ(''\--.q)u 



-t-..g(i°)[4A(3°)(^:g)i2+A(i2)(^:g)i2] 



+2(.Q(io))M2t'^'^^^'^^^^^ ■ b-(s^:9)]i2 



(ei2)^ 



(ei2)2 



}.(E14) 



If we set rii = 71,2 in eqs. (|E12IIE14| . the powers of ei2 
in the denominators will be canceled, thus rjj ci^e 

finite in the limit ei « 62, as they should be. Moreover, 
if we set 62 = ei> = ei in the resultant expressions, we 
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obtain eqs. (IE6IIE8I) . derived from the saturation princi- 
ple. This is also true in the case of a general Vl234- With 
the solutions (IE6p and (IE12I) we can calculate Q'^'^°\ 

r)(20) _ V^ „(20) V^ (20) ^ 

y Z^ ^12 121 + 2^ r-^y qvi - 



KQi^O) 



ei#e2 



ei2 






where we have used eqs. (jE4p and (IE5[) . Canceling 
q(20) fj-om both sides we obtain eq. ([S^]). Similarly from 
Q(02) ^ ^^^__,^^ ''12^^921 +Eei=ei, ^ir^9i'i ^c obtain eq. 

In the quartic order, the leading ei = ei' matrix ele- 
ment T^y is determined from eq. (P5)l . 

2zr(f - 2zA(3°)r[°?) -t- 2^A(12)^(20) = 
-«Q(20)[^^ ^(01)]^^, _ 2kQ('%, r(")]n', (E16) 



with the solution (ei = ey) 



(30) 



''ii' 



-2(«Q(io))2a(30)^ 

2 

-2(a.Q(i"))2A(12)^ 



gi2g2i''^i2 



(ei2)^ 



qi2q2vni2 



(ei2)= 



(10),.n(20) V- gl2g21'»i2 . ^(10)r rii) 



-«;Q(io)Kg(2") ^ 



(ei2)2 



+ iKQ'-'^'[q,r^''>]n'm7) 



.(30) 



The leading ei 7^ 62 matrix element r-fj is determined 
from eq. ([57]) . 

- 3zA(3«)rf,i) - 3*A(4°)r(°i) - ei^r^f ^ + ^^g(^'')ni29i2 
with the solution (ei 7^ 62) 



(11) 
"12 

ei2 2' 



_3^A(30)!i^ „ 3^Q(10)^g(20) Mil^ 



ei2 

^3^Q(ro)M!%2, 
2 ei2 



(E19) 



Then from Q(;;o) = Ee.^e. '^r'?2i + Ee,=e„ r[f!>qvi we 
obtain eq. 



The solutions r^™" are needed if we want to calculate 
the transitions of the operator alai from eq. (ITil) . 



Appendix F: Quartic Potential Dominance 

Around the critical point w^ w the stability of the 
system is restored by higher order anharmonicities. We 



assume that the quartic potential term A'^^^^a'^/4 is dom- 
inate, and study the conditions for this to be true. Under 
the rescaling of a and tt 



= (A(40))Ka, 



(A(40))- 



(Fl) 



which preserves the commutation relation [a,7r] — i, the 
Hamiltonian ^ is written as 



i7-£;o = (A(4o))K(^« 



1-4 , 1-2 

—n 
2 



A(30)(A(4"))-f 



A(12)(A(40))-i A(22)(A(40))-i 



A(04)(A(40))i_^ A(6")(A(40))-| 



.(F2) 



-a'* is dominant if coefficients of other 



gives 



Thus the term 

terms, e.g. A(^°)(A^^°^)~3 , are small. We consider their 
dependence on the collectivity factor J7 in the factoriz- 
able force model. Let the quadrupole operator q have 
the property specified in Appendix [D] Eq. (IE5|) gives 
(kQ(i°))2 ~ 0-1. Eq. dH gives A^^o) - ^-3. Eq. dH 
A(i2)-|-2^(^ ~ ri"^, and we assume A^^^^ -- f7"3, 
- rj-5. Eq. dHSl) gives A^^o) ^ f]-i. A consis- 
tent estimation gives A^^^) ^ f^-i^ A^"-*) ^ n^^. In the 
expression of A^^^^ there should be terms like (kQ^^"-')^- 
[trace with six g's], thus A'®°' ^ fi"^. In conclusion, 

A(30)(A(40))-i _ ^i^ A(i2)(A(«))-^ ^ n~^, 
A(22)(A(40))-i _ ^-f ^ A(04)(A(40))§ ^ f^-f , 

A(60)(A(40))-|_j^-§_ (F3) 



The estimates (jF3l) are consistent with those in Ref. [22| . 
All terms except A'^'^^^ are suppressed by powers of l/fi. 
The A(^°' term is given by three-body loops (|83p . which 
are usually suppressed, because of cancelations due to 
the approximate particle-hole symmetry near the Fermi 
surface, similarly to the Furry theorem of QED. In the 
case of a spherical nucleus, A^'^^^ should be small. 



Appendix G: Details of Realistic Nuclear 
Application 



Here we supply the details for Sec. |Vll In eq. ([90)) 
R is Hermitian, K is antisymmetric. The Hermitian of 
K is {K'')i2 = ala\. W{R} and f{R} in eq. ^ are 
Hermitian, A{K} in eq. (IM|) is antisymmetric. The 
Hermitian of A{K} is 



A^{K},2 = (A{i^}2l)t = iE^4312(i^t)34. 



(Gl) 



34 



The expansion of the operator R replacing eq. ([M)) is 



R = p + R 



N 



E^^" 
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L=0,2,4 pi L=0a,2,3,4 M 

4 E EC(-^xA' +^ E E E^SW{("^x"^)''"%" 

L=0,2,4 M L=0,2,3,4,6;=0,2,4 M 

4 E ECW,-^}.^ +J E E E-SU{(-^x.ty,.t}. 

L=0a,2,3,4 M L=0,2,3,4,6i=0,2,4 ^t 



^ E E'ri("^x"'^'"^}" +iE E E-SU{(«^xat)vt}^ 

L=0,2,3,4,6 M L l=0,2A fJ. 

I E E^S^{(-^x.ty^.t}^. +1^ ^ E-iS{i«^-^}''"^}.' 



L=0,2,3,4,6 fJ. L 1=0,1,2,3,4 fi. 

+i E E E-S{("^xaty,.t}. +1^ ^ E«i234W,(-^x.ty}. 

L=0,l,2,3,4,5,6 ;=0,2,4 M i ;=0,2,4 M 

+i E E E-Sw,(-^ x.ty },- + ... (G2) +i^ 5: E''Sl°3l«"^-^}^-^' 

L=0,l,2,3,4,5,6 (=0,2,4 M L /=0, 1,2, 3,4 fi 

Three identical d bosons can couple to L = 0, 2,3,4, 6. +■•■ l"-^Oj 

In the a3 and n^ terms of eq m we choose the in- j^ ^^_ ^^^ ^^^ ^3^ ^2^^ „^2 ^^^ tt^ terms are over- 

termediate quantum number for each L to be h; this complete. This form is convenient for finding expressions 

choice does not mnuence the results, it is Hermitian, („„) . („„) , , . • • i 

time-even, invariant under rotation and parity Il2|. This ^'^^^ . , . . \^ rm o- -i i i i 

(•„„) as explained m Appendix |H| Similarly we need the ex- 
implies that the coefficient rV,, has the same symme- . r atf t i mt i atf t t t i j 

^ , s pansions oi iV [a4a3a2aij, i\ [a4a3a2ai\, l\ [a!^a!^a2ai\ and 

tries as the operator part {a™, tt"}: (r)™"^)'!' has angular Nla^a^a^a^]. 

momentum L and projection fi, even parity, sign of (— )" 

under time-reversal, {rj^" Y — (~)^^^^i™u • Similarly 

the expansion of the operator isT is 1- Exact Equations of Motion 

K = K + K^ = K + Y^ kl^^'fal + ^ fc^^^^t The Hamiltonian (dJ in the normal ordering form is 

+ - E EC("'xat)^ + ... (G3) H^mH\<^)+Y^f{p},,N[aU 



L=0,2,4 M 1 1 

+ - "V Al'{K}i2A^[aia2] -I- - "V A{K}i2A^[a|4l 
K is anti-symmetric, time-even, invariant under rotation 2 -^^ 2 ^-^ 

and parity. Thus (fc^^" )^ has angular momentum L 2. 

and projection /i, even parity, sign of (— )" under time- +7 2^ ''1234 jV[a]^ 020304], (G7) 

(mn) _ _. (mn) rp, „ U„^^; + ;a„ r^f ^r. irrr^ 1234 



reversal, fcl 12 ~ ~^Lu2i- ^^'^ Hermitian of eq. (jG3 

is where 



^' = ^' + (^^)" - ^' + E(^^)r "1. (<i>iiii<j>) = 5:(z,2 + V{p}i2)P2i 

^ 12 

-TjJ^)r< + \ E E(^^)ff(-^x«t). + ...(G4) +iEAWi2(-t),^ (G8) 



M L=0,2,4 /J 

where jg ^^j^g average energy on 1$). The exact e.o.m. in the full 

fci^ ^ (-)^-^fcL-p ^ fcL = (-)^"^4-^- (G5) 'P^'^^ replacing eq. m is 



''L/i V 7 L — fi 

The expansion of A^[a4a3a2ai] replacing eq. (jisp is 



t„t„ „ 1 !_;„„_ /rra ;„ [i?i2,ii] = [alai,H] 



[f{p},p]i2 - {^ At{K})i2 + (A{4 ^«t)i2 
7V[al4a2ai] = i ^ E^S234(«^ x at)L +[f{p},R%2 + [W{R''}, p]^2 - {K'' AHk}),^ 

i=0,2,4 M -(k At{i^^})i2 + {A{k} iftW)i2 + (A{A'^} «t)^2 

+ 2 E E^Si234(^^x^t)i +lY.^V,345N[4alaia5]-V5i32N[alala3ai]),iG9) 

L=0,2,4 M "^ 345 
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and 

[K^2,H] = [a2ai,H] = {k f{p))i2 + (/{p} '«)i2 
+A{4i2 - (A{k} p^)i2 - {p A{k})i2 

+ (/^^ /^{p})l2 + {f{p} if^)l2 - (A{4 {R'f)l2 

-{R^ A{4)i2 + (a. W^{R^})i2 + {W{R^} k),2 

+A{if^}l2 - (p A{if^})i2 - (A{i^^} p^)i2 

+ - ^(V'2543^[a5°^4a3ai] - V^i543^[a5a4a3a2]).(G10) 

345 



Then in eqs. (|G9p and (JGIOP we equate the l.h.s. and 
r.h.s. coefBcicnts of the same phonon structure: 1, a^, 
TT^, ^{a X g)^ . . . , and obtain e.o.m. in the collective 
band (l97lfT06ll of Sec. lyTBl 



2. Hartree-Fock-Bogoliubov Equation 

The HFB equation ^ says that S^™) and £)(°°) can 
be diagonahzed simuhaneously, 



E 
-E 



n 
/-n 



0, (Gil) 



where E and n are diagonal matrices. The chemical po- 
tential p (buried in /) is determined by iV = J2i Pn — 
Tr{p}. The unitary canonical transformation from the 
original s.p. operators a[, oi to the new quasiparticle 
operators b\, b\ are 

bx = ^{ulxai - vixa\), b{ = ^(uiAflj - i;i^ai).(G12) 

1 1 

If |<i>) is a "quasiparticle determinant", 1$) = norm • 
J|^6a|0), then the normal ordering with respect to |<I>) 
is to put fot's to the left of 6's. Eq. ^ gives: 



-vu 



(G13) 



In this case Z3('^") is diagonalized by the canonical trans- 
formation (|G12p : 



C/ = 



-V' u 



UD^'^^^U^ 




/ 



(G14) 



where the matrix n in eq. (IG11[) vanishes. The HFB 
equation ([W)) requires 5^°"^ is diagonalized by U simul- 
taneously. In this article we assume that |<i>) is a "quasi- 
particle determinant" . 

It is convenient to solve the e.o.m. (19711106^ in the 
quasiparticle basis (multiplying U from left and t/^ from 
right). The density matrix operators in this basis are 

Rl^ = blbi = N[blbi], Kl^ = b2bi = N[b2bi], 

{K'^h2 = blb\^N[blb\].{G15) 

R^ is a mix of R, K, K^ of eqs. ^■, so do K^ and 
K^^ . The expansions of them are defined similarly to 



eqs. (IG2IIG4P : 

The field matrices in the quasiparticle basis are 

D^ = UDU\ S^ ^ USUI (G17) 

We need to express them in terms of i?'', K'' and X^^ of 
eq. ((Gl6|) . The result of D^ is simple: 



D^ = UDU^ = 



R" 



X" 



(G18) 



Sa Sb 
Sc Sd 



(G19) 



The result of S^ is long: 
S^ = USU^ - 

where 

SA = u^f{R}u-u''A{K}v* 
-v^A'fiKju - v^,f{R}v* 

= ( u^Zu - v^Z^v* + v)W{vv^u 
+ ( u'^W{u.p{u'^)a.}u + u'^A{u.ii[v'^)a.}v* 

-v)W{v.a{v'^)p.}u - v)/\{v.a{u^)l}.}v* 
+ i - v)W{v.p{v))a.}u - v)A{v.fi{v^)a.}v* 

and 



Sb = -u\f{R}v + v)A{K}u* 
WA'^{K}v + v^f{R}u* 

= ( -u'^Zv + v^Z'^u* -u'^W{vv^}v 

-u^ A{vu^}u* - v'^A'^{vu^}v -f v'^W^{vv^}u* ] 

+ f - u'^W{u.p{u'^)a.}v - U^A{u.;3(f^)a.}u* 

+v)W{v.a{v'^)p.}v + v)A{v.a{u^)l3-}u* 
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+ ( u''W{v.f){u'^)a.}v + u'^A{v.piv'^)a.}u* 

+ ( u''W{u.f:!{v'')a.}v + u''A{u.f){u^)a-}u* 
+v^A^v.^{v^)^.}v - v'^W^{u.p{v^)^.}u* )6„6^,(G21) 

and 

Sc = --v^f{R}u + v'fAiKjv* 
+u^A^{K}u + u^f{R}v* = {Sb)\ (G22) 

and 

SD=v''f{R}v-v'^A{K}u* 
-u^A''{K}v - u'^f{R}u* = -(5"^)^. (G23) 

From now on we will always work in the quasiparticle 
basis unless otherwise specified. For simplicity we will 
drop the superscripts '' in eq. (jG16p and ^ in eq. (|G17p . 



3. Quasi-particle Random Phase Approximation 

The QRPA equations ^E^ and ([Ml) are (we have 
dropped the superscript ^) 

Each of the above two equations has four components, 
only two of them are independent. The upper-left com- 
ponent gives 



(10) _ (01) _ „ 
Vl2 ^ Vl2 ~ ^' 



The upper-right component gives 



*'*^/j12 — WBV12 



-iu! k 



2^,(01) 



/il2 



(Sb) 



(10) 
/il2 



(Ei- 
-{Ei 



E2)k^, 



The formal solution is 



^coHSBC^~{E^+E2){SBi\'^ 
(El + E2y - uj^ 



^(01) _ -(-El + E2){Sb)]^12 ~ iiSB)f,i2 



^m) 



(10) 



>12 



[El + E2)^ - uj^ 



(G24) 

(G25) 
(G26) 

(G27) 
(G28) 



From eqs. (IG2T1) . ((G27t and (1628]) we obtain a linear 
homogenous set of equations for (5'b)^^*'' and {SbY°^\ 
a non-zero solution requires a zero determinant, from 
which we solve for o;^. 



Again to fix the normalization of k 



(10/01) 
>12 



we need the 



saturation principle. Since now we are solving everything 
in the quasiparticle basis, it is convenient to redo the 
saturation principle in the quasiparticle basis. After that 



we obtain the normalization condition (independent of 



■ r t ] constant terms 



I Tr{(fc(°i))t kl^^n + I Tr{fc("i) fctdO}. (g29) 



4. Cubic Anharmonicity and Quartic 
Anharmonicity 

The second order e.o.m. are eqs. (|100m02|) . Dj^Ji 3.„ 

is determined from eq. (|102p alone. Dj^^q 4."^ is deter- 
mined from eqs. (|100m02|) . They are expressed in terms 
of lower order quantities. When L — 2, A^'^") and A^^^^ 



,(20/11/02) 



is determined in 



enter eqs. (|100m02l) . and D)~^'^.^ 
termsof A(30) and A(i2). 

Similarly to the situation in Sec. IIIDl the Ei = Ey 
matrix elements Tj^Jq 2 4 uii' ^^^ ''l=i 3uii' ^^^ fixed by 
the saturation principle {Ei — Ey ) : 

L = 0,2A: ra,=2(fc(io)xfct(io))L^^,, (g30) 
L = 0,2,4: rf;],, = 2(fc(oi) X P(°i))^,„ (G31) 

L = 0,1,2,3A: r[%,^ik(^'^xm^Xn' 

+(-)^(fc(oi) X (fct)(i"))^^^,. (G32) 



Eqs. (|G30IIG32t are consistent with the second order 
e.o.m. (|100m02p . 

The third order e.o.m. are eqs. (|103m06|) . The L ^ 2 
quantities (-D/S*)^ / J are solved in terms of lower 

order quantities. When L = 2, A^^o), Af^^ and A^o") 
enter into the equations and we have the solvability con- 
dition as explained in eqs. (|109p and (|110p . 



Appendix H: Saturation Principle for Section IVII 

Keeping only one-body terms in a^ and tt^, 

a^L = ^(a:Mi2ai"2 + 2^*1201^2 + zli2^{al), (HI) 
12 

^!^^L = ^(PA'i2aia2 + o^i2aia2 + oli2<^{al), (H2) 
12 

where 

xl = i-Tx^,, pI - (-Tp-,, 

^/il2 = —Z^21, 0^12 = —0^21, (H3) 

we have the following identities in the full space: 

[i?i2, a^] = [x^,, R]i2 - 2{Kz^)i2 + 2(ztxt)i2, (H4) 
[Ki2,a^] = {x^K)i2 - {x^K)'^^ 
+2zl,, + 2{Rzl)'^,-2iRzl)i2, (H5) 

[iK^)i2,a^] = -{K^x^)i2 + iK^x^)^2 

-2z^i2-2(z^i?)f2 + 2(z^i?)i2, (H6) 
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and 

[i?12, TT^] = [p^,R]i2 - 2iK0^)l2 + 2{olK^)i2, (H7) 

+2ol,2 + 2{Rdl)'^2-2{Rol)i2, (H8) 

[iK^)i2,7r^] = -{K^P^)i2 + {K^P^.)l2 

-2o^i2 - 2(o^i?)[2 + 2(o^i?)i2. (H9) 

We obtain a set of equations by equating the l.h.s. and 
r.h.s. coefficients of the same phonon structure: 1, a^, 
T^fj., ^{a X a)j^ . . . Considering the length we do not hst 
them here. 

Similarly we calculate the commutators of 
N[a\ala2ai], N[a\a:ia2ai] and A^[a4a3a2ai] with 
a^ and tt^. We give only the result of [Af[a4a3a2ai],a^] 
as an example: 

[N[alala2ai],af_,] = 

+^r%'^N[a\a2] ~ ^k%'^ N[a2a,] - ^kl^['^ N[alal] 

- ^ Xf^25N[alalaia5] + ^ a;^i5iV[aJ4a2a5] 

5 5 

- X! Xt^54N[alala2ai] + ^ Xf^53N[alala2ai] 

5 5 

+2 ^ Zf^54N[ala2aia5] + 2 ^ z^35^[a4a2aia5] 

5 5 

+2 E i5iN[444<^2] + 2 J2 ^i25^[444«i],(Hio) 

5 5 

where we have used the lowest order results from eqs. 
(|H4IIH9p . Equating the l.h.s. and r.h.s. coefficients of 
the same phonon structure we obtain a set of equations. 
We give only the (a x a)^ terms as an example. Using 
results from eqs. (JH4IIH9J1 we have 



4 E E 



V2i + 1 



^V5 



L/=0,2,4 L 



H'l)i234 X (a X a)^f) 



+ (^i23 X '^U ) + (^(12 X fc34 ) + (r;^4 X r23 ) 



(-.^3^ X ^^ro" + (fc;3^4^^ X k\^'>r i x (^ x a)^ 



(Hll) 



The l.h.s. and r.h.s. of eq. (|Hlip come from the l.h.s. 
and r.h.s. of eq. (|H10I) . respectivelv. The following ex- 
pressions satisfy eq. (|Hlip : 



(21) |L=0,1,2,3.4,5,6 
'tn234 U=0.2,4 



(20)^(01)nL 



(^il4 ' 23 
(20)^(01)^L 



r + (r. 



(20)^(01)xL 



ir. 



(20)^(01)^L 



(13 ' 24 



-,(20) 71(01) xL 



:t(20),(01).,L 



V'i24 'l3 ) '^ \'^112 '^34 ) '^['^134: "'12 



r,(Hi2) 



(11x10) |L=0,1,2,3,4,5,6 _ 
^L11234 U=0,l,2,3,4 ~ 



2-{C 



.(11) ^ ^(lo)^i 



(14 ^ '23 



, /„(") y „(10)^L / (11) y ^(10)nL / (11) ^ „(10)nL 
"•"^(23 ^ '14 ) \'ll3 ^ '24 J Vi24 ^ ' 13 J 

+ (fcL"^xfcr^)- + (fctai)xfc(f)M.(H13) 

Appendix I: Values of gf'i, and ■yf'i, 
The definition of gf^^, is given by eq. (llOSp . 
i{{a,7r}'',a}^= E 5m' ' ^{(« x a)',^}^^, 

;=0,2,4 

which implies 

5,^, = (-)'-V(2; + l)(2;' + l)| 2 J ^/ } =5^^(^) 
The definition of 7/";/ is given by eq. (jl07p , 

{(a X a)', a},^ = 7;^;' ' {(" x a)'', a}^. 

Analytical expressions of 7;^;, can be obtained in the fol- 
lowing way. Assume / and I' are even. We have the 
identity 



vt V r.t\l' ^UL 



1 |L=0,2,3,4,6 



{{a^xa^y^aX. ('t x 7r)!J ] |--' 



2,4 



':r7F E E E(-)'^^' \/(2i + 1)(2; + l)(2L^- + 1) ^^ E ('^M'+2-g,^. )-{(«^xat)',7rt}^. (12) 



Lf=0,2.4 L 



2 2 / 

2 i i/ 



('■m234^ X (a X a)^/) 



i E E 



V2L + 1 



V5 



Lf=0,2,4 L 



([-(<iVr) 



+(^L^i4''23 '» +(ryA^r 



^^ , /^(20) „(01)xL _ / (20) (Ol)xL 
/ ^"V'L/23'l4 ) \'Lf24'l3 ) 



i=0,2,4 

Replacing V in eq. p2)) by I" we obtain 

[{(atxar^,t}L^(,X,)0]|i;=J)^2^.3.C6^ 

*4 E ('^M"+2-.g,V)-{("^xat/,vrt}^^. (13) 



(20) (Ol)xL 
L/14'23 I 

+(4%fcr)'+(^i%"Mf)"]x(«x. 

-*7f E E E(-)'^^' V(2L + 1)(2; + 1)(2L/ + 1) 



'V5 



1=0,2,4 



Lt=Q,2.4 L 



2 2 I 
2 L Li 



}-([-(4^^^x.if)- 



Let us take the ratio of eq. (|l2|) /ea. (|I3| . The l.h.s. is 
^f!i„ by eq. (fT07| . Since in the r.h.s. of eq. dH]) and 
eq. (jlS]). {(a^ x at)',7r^^}^ with different ^ are linearly 
independent, we have 

(S; /' + 2 • 0/',, 
^'V = r-— y4^' G = 0,2,4). (14) 

d/,r' + 2 • 5/;;,, 
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The ratio on the r.h.s. is independent of /. Since the 
matrix Su' + 2 ■ gf^, is symmetric (with respect to I, I'), 
eq. p4)) imphes 



(15) 
where ff^ = Jl + 2 ■ gf^. Then from eq. (HH we obtain 

lv,v' = tI- (16) 

We will use only L = 2: 



fL=2 _ / ' fL=2 _ ^ fL=2 _ " (T^\ 

In the main text the superscript ^^^ on //"^^ is dropped 
for simplicity. 



Appendix J: Conventions 

Our convention for the Wigner-Eckart theorem is: 

{n^3im,\T^\n232m2) = C]lZlx^. ' ("iJi||r^ll"2j2). (Jl) 
The reduced matrix element of q\ (|112p is 
WWllWnl'j') = 

d' + \-l(_\i' + \-i / (2A+l)(2j' + iy 



47r(2j + l) 



■Vi,;,o/'^^'^'^^(^)^"'j('^)^"'''^'('')' ^' + A-nseven, 
0, /' + A - Hs odd, 

where the s.p. levels \nljm) are defined as 



YTiljm ^nlj V j 



C^Z\ „™ i'Yi^i[U,(i)) Xm,, 



/ ^ Imi.sm 
rainis 

in which spin s = 5, Rnij{r) is a real function, and a 
factor ?' is included. 

In this article we have used the matrix elements of the 
realistic quadrupole moment operator in the harmonic os- 
cillator s.p. basis. In this case A = 2, f\{r) ~ r^, Rnij{r) 
is independent of j, and n = 2nj, + I is the major-shell 
quantum number. The non- vanishing matrix elements of 
(nljWqmn'l'f) have n-n' = -2, 0, 2, and l~l' = -2, 0, 2. 
For these combinations the symmetric radial integral be- 
comes 



dr r'^Rni{r)Rn'i'{r) 



n + 



n' = n, V = I, 
n' ^n, I' ^1 + 2, 



h'-l 



~\^ J {n + I + i){n - I + 2) , n' = n + 2, I' = I, 
^y /{n + l + 3){n + l + 5) , n' = n + 2, I' = 1 + 2, 
, ^y/{n-l + 2){n-l + 4:), n' = n + 2, I' = I ^ 2, 



where b = \/—jr- is the length parameter, fig is the har- 



mOo 



monic oscillator frequency. As mentioned at the begin- 
ning of Sec. I VII B[ the factor b^ will be combined with k 
to make qj^ dimensionless. 



Appendix K: Details of Quadrupole plus Pairing 
Model 



Here we supply the details for Sec. I VIII In the 
quadrupole plus pairing model the HFB equation be- 
comes the BCS equation. The canonical transformation 
(|G12|) becomes 

Wl2=<5l2Ml, W12 = -(Jial'l, (ui)^ + (wi)^ = 1,(K1) 

where ui = uj, wi = fj are real numbers. The density 
matrices (|G13p become 



P12 = <5i2(i'i)^, K12 = Sj^^'^ivi. 
The /{p}i2 field dH]) becomes 

f{p}i2^Si2[ei-fi-G{vi)^] 



(K2) 



3 M 



61261 - K(gi^o) 



12< 



)(00) 



(K3) 



where Q^™) = Y.3ili^=o)33{v3)^ , and ei = ei - /i - 
G{vi)'^ = ej, and we neglect the incoherent sum. In 
the case of a spherical mean field, Q^^^^ = 0, thus only 
the ^1261 term survives. The A{k}i2 field ([^^ becomes 

G 

A{k}i2 = -17 ^uYl "3^3 ~ '*X1 XI ^Il3'?^23 "3^3 



3 /J 



-S,^A, (K4) 



where the pairing energy A = ^ S3 ^3''^3' ^^'^ ^^ neglect 
the quadrupole-force contribution to the pairing poten- 
tial. The HFB equation (|97| gives the BCS set of equa- 
tions pi6m2op . 

Quadrupole moment in the quasiparticle basis is given 
by 

12 

^(wi)^g^ii -f- ^(uiU2 - viV2)qf_i2iblbi 
1 12 

--^{uiV2 + U2Vi)qf,2i{b2^i + blbl). (K5) 



12 



Substituting the expansions of i?^ K'' and iC^t ((GT6)) 
into eq. (|K5p we have 



^^,3„, {(axa)'(^=2),a}^ 



6 



■E • 

Z=0,2.4 



^(i2){a,(7rX7r)'}^ 



(K6) 



where Q^™") is expressed in terms of r's and fc's in eq. 
(|G16p . Note on the r.h.s. only terms with the same 
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symmetry as Q^ survive. The Hermitian property (J113p 
implies that ah Q^™") are reaL Pairing operator is given 

by 

1 1 

+ ^ uiz;i6l6i + \ Y}^u,)\b^ ~ (z;i)26l5t]. (K7) 



[(£;^+i52)2_^2]2 «||12^||2r l^i4J 



The cubic potential term (|124p in the original form is 

A(3°)=-(.g(i°))MTr[(C«xe(^))^=\tj 

+Tr{77,f(i)}^=2CWtl}^ (K15) 



Substituting the expansions of i?^ iiT'' and K'''^ (|GT6l) 
into eq. (jK7l) we have 

(a X q:)o 



P=-i^7.iz;i+p(20) 



2V 2 



p(ii) {«:7r}g „(02) (tt X 7r)g 
2 2 



where each term on the r.h.s. is real and independent of 
^. 

The quartic potential term (|126p can be written in the 
form of reduced matrix elements. The result is long, we 
write down them term by term. The term with (A^'^"')^ 
is 



(K8) 



P + pt is Hermitian and time-even, P — P^ is anti- 
Hermitian and time-odd. Thus p(^°) and p(°^) are real, 
P^^^^ is pure imaginary. The f{R} field (|9ip becomes 

f{R}i2 = <5i2ei - GN[ala~^] -^Y^ql^^Q^^ 

34 p p 

Here again we neglect the incoherent sum and the pairing 
contribution beyond the mean field. The pairing field 
becomes 



-2/2-(^Q(i°))2(A(30))2Tr 



£ f(5)t 
(A(30)) 



The terms with A(3°)(kQ(i°))3 are 



/2-A(30)(«;Q(io))3( Tr 



-Tr 



j(3)i.L=2<r(l)t 



{^/,e(^^}rei 



f 2 Tr 
-Tr 



"/„ ;^(l)li=2^(3)t 



A{i^}i2 = <5i2GP - A. E E 9ii4'?p23a3«4 « 5,^GP, (KIO) 

34 p 

again neglecting the quadrupole-force contribution. Fi- 
nally, the field Sa/b/c/d (!G20IIG23|) become 

{SA)i2^5^2[{uif-{vif]ei 
-n[uiU2 - W1W2] E 9ii2<3M - uiVi6i2G{P + Pt),(Kll) 

(S'b)i2 = 2wiwi(5i2ei - n[uiV2 -f U2«i] E 'ilu^f^ 

+<5i2G[K)2p-(t;i)2pt],(K12) 



f(2)j.i=2^(2)t 

= 2A(30)/2.(,.g(i"))3 J2 



nijl"2j2n3j3 

(2ji + l)(2J2 + l)(2j3 + l) f 2 2 2 

5 \ Jl J2 J3 

4f(l)M3)t t , f(2)t^(2)t t 

^||13 ^||32 ''||21 ^ ^||13 ^||32 ''||21 



(K17) 



and Sc = {Sb)\ Sd = -{SaV- Substituting eqs. (lK6l 
and (jKSp into the above equations we obtain the expan- 
sions of Sa/b/c/d- ^^ 

The QRPA secular equation (|122p in the form of re- 
duced matrix elements is 



The two terms with p(20) are 

- /O • GP(2°)(«Q(1°))2 E("l^l + "2^2)^.12^1^2^ 
12 

= -/o-GP(20)(«:Q(io))2 ^ (- 



«2-Jl 



nijin2j2 



V(2ji + l)(2j2 + l) 



■ {uivi + U2V2) ■ ^l^^l^W {K18) 



and 



1 = , ^ ( )..-.. V(2j-i + l)(2j2 + l) 



"ljin2j2 



-£'1 + -^2 



(^1+ ^2)2-^^2 '-||12'-| 



ai2a2r (K13) 



The normalization condition (jl23p in the form of reduced 
matrix elements is 



i = (.g(i")f E (- 



'\32-3 



^V(2ji + l)(2j2 + l) 



"Ul"2j2 



+ /o-GP(20)(Kg(io))2 

v-r M^^M^ , M^_M^. ,(1) t 

•Z^i 2£;i "^ 2^2 ^V12^m21 

12 ^ "^ 

= /o-Gp(2")(Kg(io))2 

^ ( ^,,„,, V(2ji-H)(2j2 + l) 

, K)^-K)^ , M^-(^2)^ , ,t .(i)t 



2Ei 



2E2 



}-'?ii2q2^(Ki9) 
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where p(^") is given in eq. (|125l) . The terms with 



1=0,2,4 



-Tr 



{^,{r;,e(^n'^(^)}rc«t 



-Tr 



[{^(l),{,,^(l)}ai)}L=2^ 



/=0,2,4 nijin2J2n3J3n4J4 

{2ji + l)(2j2 + l)(2j3 + l)(2j4 + 1)(2; + 1) 



''||13?||34'5||42 ?||21 ^'l\\ 



5 
/ 2 2 

jl J2 J3 

1 



113 E.+E, 'II34HI 



2 2 ; 

J3 ^2 J4 



„t ^(l)t^(l)t 
''||34S||42 ?||21 



-2^L3 • 



1 



f(l)t^t t(l)t 



113 K. + i;. 



ei;;i^S42Cir2;' ■(k20) 



Q ■ Q part, k is the strength of the octupole force. The 
mean field is determined by the HFB equation. In the 
harmonic order the two modes do not mix, the octupole 
mode satisfies the same QRPA equation (|122l) and nor- 
malization condition (jl23p as the quadrupole mode, with 
necessary changes. In the next order we have the main 
result: 



(10)n,2 



E 



A(i°I^°)=2kQ(i")(«Q 



(2ji + l)(2j2 + l)(2j3 + l) 



ni]in2j2n3]3 



7 



3 2 3 

ji J2 id. 



\ . r 2 f (i)t £(i)t ^t , t p(i)t ^"(i)t 

r [ ^||13 ?||32 ''||21 ^ ''||13 ?||32 '5||21 



.(L3) 



The octupole operator q connects s.p. levels with oppo- 
site parity, thus the intruder state becomes important. 
This may destroy in eq. (jL3P symmetry with respect to 
the Fermi surface. Three-body forces will contribute to 
the A(i°I^°) term quite differently. 



Appendix L: Mode Coupling 

In many soft nuclei there exists a low-lying octupole 
(3^) mode. It can interact strongly with the quadrupole 
(2+) mode, and both of them should be kept in the col- 
lective subspace. For convenience we still use a^, tt^ for 
the quadrupole mode; and use d^, tt^ for the octupole 
mode. The collective bosonic Hamiltonian replacing eq. 
© is 



H = Y\/5(a X a)l + ]^Vl{^ x ^)l 



U! 



H — —V7{a X a 



lg + iV7(^x^)° 



A(10|20) 



V7{a X {ax a)^)o + 



(LI) 



jY(io120) jg ^jj^g most important mode-coupling term in the 
case of soft vibrations with large amplitudes. Following 
the procedure of Sec. |ll] and IIIIl we are able to deter- 
mine the leading constant term of A'^^'^l^*'^ in a Taylor 
expansion over both w^ and w^ [see eq. (|47|]. Below we 
give the result in the quadrupole plus pairing model. The 
microscopic Hamiltonian is: 



G 



H = E(ei - ^^■)a\al - J E aiai«2a2 



12 



1234 /I 



(-'«'Zu14<?m23 + l^li 



13'?A'24 



-«9ii49M23 + i^ql^:i%2'i)a\a\a:>,ai. (L2) 

Approximately, this Hamiltonian can be written as if ~ 

the difference is in a one-body term originating from the 
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TABLE I: The first excitation energy Ei — Eq for different J in tlie Lipkin model. The first three lines are the results of 
diagonalizing different bosonic Hamiltonians in different phonon spaces; the last two lines are the results of diagonalizing eq. 
(|58|l directly in the {\JM)} space. In the space \n < 2J), the matrix of e.g. a* is calculated by multiplications of the a 
matrices, which is different from truncating the a"^ matrix of the space \n < +cx3). Higher excited states from the GDM method 
are also in good agreement with the exact results; please see the figure in Ref. [23 |. 



El ~Eo 


J = 2 


J = 3 


J = 4 


J = Q 


J= 10 


J = 50 


Isl + ll - J-Hl m Ir? < 2 n 
J 4 + 2 8.J 4 ™ |n ^ zj; 

5^ + ^in|n<2J) 


1.011 
1.087 
1.087 


0.904 
0.949 
0.950 


0.835 
0.863 
0.863 


0.740 
0.754 
0.754 


0.630 
0.636 
0.636 


0.371 
0.372 
0.372 


exact , K = 2j+T 


0.950 
0.895 


0.869 
0.776 


0.808 
0.707 


0.722 
0.625 


0.620 
0.537 


0.370 
0.334 


l^j'vLikjL ^ 1 'J or 1 
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TABLE II: Results of the quadrupole plus pairing model at different pairing strength G. All quantities are in unit of MeV. A 
is the solution of eq. (|116p . The chemical potential /i is the solution of eq. (|120p . k.c is the critical k such that uP in eq. H122p 
becomes zero. A"''' is given by eq. ([124]). A'*"' is given by eq. H126P setting A'^^' = 0. "GDM E2+" is the excitation energy of 



the first 2^ state by diagonalizing eq. (|96|) for : 



A(12) ^ ^(04) ^ ^ 



(22) 



0. "NuShellX E^+" is the exact excitation energy 



of the first 2'*' state by diagonalizing eq. IjllSp . in which G and «; are given by G and Kc in the table. Similarly "NuShellX 
-E^+" is the exact excitation energy of the first 4^ state. 



G 


0.03 0.06 0.09 0.11 0.12 0.15 0.18 0.21 0.25 0.30 


A 


0.0 0.0 0.0 0.0 0.0 0.066 0.453 0.672 0.862 1.096 1.374 


M 


0.5 0.5 0.5 0.5 0.5 0.454 0.444 0.429 0.415 0.395 0.370 


K.C 


0.102 0.105 0.107 0.110 0.112 0.113 0.113 0.122 0.135 0.154 0.179 


A(30) 


-0.160 -0.173 -0.188 -0.203 -0.213 -0.219 -0.234 -0.270 -0.310 -0.378 -0.474 


A(40) 


0.483 0.526 0.572 0.621 0.655 0.616 1.185 1.918 2.901 4.683 7.830 


GDM £2+ 


0.882 0.908 0.933 0.959 0.976 0.955 1.194 1.405 1.614 1.894 2.249 


NuShellX E^+ 


0.855 0.892 0.944 1.023 1.106 1.158 1.353 1.552 1.764 2.059 2.438 


NuShellX E,+ 


0.778 0.827 0.927 1.110 1.284 1.383 1.705 2.076 2.465 2.987 3.631 
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FIG. 1: (Color online) The first excitation energy Ei — Eo in the factorizable force model, as a function of k. The red circles 

- in the infinite phonon space. The black triangles give the RPA frequency lu, this 



u j a 
2 



+ 



result from diagonalizing 

corresponds to diagonalizing ^^^ — h ^. The black squares, green crosses and purple pluses are the exact shell model results 
with different values of a and b (two lines "a = b = 0" and "a = b = 0.353" closely overlap and are indistinguishable on the 
figure). The second excitation energy E2 — Eo at k = Kc = 1/16 with a = 6 = is also in good agreement; the exact one is 
1.236, while the GDM gives 1.269. 
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FIG. 2: (Color online) The first excitation energy Ei — Eq in the factorizable force model, as a function of the parameter 
a — —b, at the critical point oj = («: = Kc = 1/16). The black squares show the exact shell model results. The red circles 

are obtained by diagonalizing ■^ + \ ■ ~ 7 in the infinite phonon space. The blue triangles and the green inverted 

triangles are obtained by diagonalizing the same Hamiltonian in two different finite phonon spaces, specified by nmax and u; u 
is the canonical transformation parameter defined in eq. (|60|) . Wmax is the maximal number of phonons. 
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FIG. 3: (Color online) Excitation energies (from Table [IT]) in the quadrupole plus pairing model as a function of the pairing 
strength G, at the critical point a; = (k = Kc). The black squares and red circles show the exact excitation energy of the first 
2+ and 4+ state, respectively, "NuShellX E^^" and "NuShellX E^+" in Table |TT] The blue triangles give "GDM £3+" from 
Table nil 
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